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Abstract. In many Direct and Inverse Scattering problems one has to use a 
parameter-fitting procedure, because analytical inversion procedures are often 
not available. In this paper a variety of such methods is presented with a 
discussion of theoretical and computational issues. 

The problem of finding small subsurface inclusions from surface scattering 
data is stated and investigated. This Inverse Scattering problem is reduced to 
an optimization problem, and solved by the Hybrid Stochastic-Deterministic 
minimization algorithm. A similar approach is used to determine layers in a 
particle from the scattering data. 

The Inverse potential scattering problem is described and its solution 
based on a parameter fitting procedure is presented for the case of spheri- 
cally symmetric potentials and fixed-energy phase shifts as the scattering data. 
The central feature of the minimization algorithm here is the Stability Index 
Method. This general approach estimates the size of the minimizing sets, and 
gives a practically useful stopping criterion for global minimization algorithms. 

The 3D inverse scattering problem with fixed-energy data is discussed. 
Its solution by the Ramm's method is described. The cases of exact and 
noisy discrete data are considered. Error estimates for the inversion algorithm 
are given in both cases of exact and noisy data. Comparison of the Ramm's 
inversion method with the inversion based on the Dirichlet-to-Neumann map 
is given and it is shown that there are many more numerical difficulties in the 
latter method than in the Ramm's method. 

An Obstacle Direct Scattering problem is treated by a novel Modified 
Rayleigh Conjecture (MRC) method. MRC's performance is compared favor- 
ably to the well known Boundary Integral Equation Method, based on the 
properties of the single and double-layer potentials. A special minimization 
procedure allows one to inexpensively compute scattered fields for 2D and 3D 
obstacles having smooth as well as nonsmooth surfaces. 

A new Support Function Method (SFM) is used for Inverse Obstacle 
Scattering problems. The SFM can work with limited data. It can also be 
used for Inverse scattering problems with unknown scattering conditions on its 
boundary (e.g. soft, or hard scattering). Another method for Inverse scattering 
problems, the Linear Sampling Method (LSM), is analyzed. Theoretical and 
computational difficulties in using this method are pointed out. 
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1. Introduction. 

Suppose that an acoustic or electromagnetic waye encounters an inhomogeneity 
and, as a consequence, gets scattered. The problem of finding the scattered wave 
assuming the knowledge of the inhomogeneity (penetrable or not) is the Direct 
Scattering problem. An impenetrable inhomogeneity is also called an obstacle. On 
the other hand, if the scattered waye is known at some points outside an inhomo- 
geneity, then we are faced with the Inverse Scattering problem, the goal of which 
is to identify this inhomogeneity, see 50, 55, 56, 17, 16) 

Among a variety of methods available to handle such problems few provide a 
mathematically justified algorithm. In many cases one has to use a parameter- 
fitting procedure, especially for inverse scattering problems, because the analytical 
inversion procedures are often not available. An important part of such a procedure 
is an efficient global optimization method, see j24L I25L I35L I36L I45L I82| . 

The general scheme for parameter-fitting procedures is simple: one has a rela- 
tion B{q) — A, where B is some operator, q is an unknown function, and A is the 
data. In inverse scattering problems q is an unknown potential, and A is the known 
scattering amplitude. If q is sought in a finite-parametric family of functions, then 
q = q{x,p), where p — (pi,....,p„) is a parameter. The parameter is found by 
solving a global minimization problem: ^[B{q{x,p)) — A] = min, where $ is some 
positive functional, and q & Q, where Q is an admissible set of q. In practice the 
above problem often has many local minimizers, and the global minimizer is not 
necessarily unique. In |55l and j57j some functionals are constructed which have 
unique global minimizer, namely, the solution to inverse scattering problem, and 
the global minimum is zero. 

Moreover, as a rule, the data A is known with some error. Thus As is known, 
such that ||yl — Aall < 5. There are no stability estimates which would show how 
the global minimizer q{x,Popt) is perturbed when the data A are replaced by the 
perturbed data As- In fact, one can easily construct examples showing that there 
is no stability of the global minimizer with respect to small errors in the data, in 
general. 

For these reasons there is no guarantee that the parameter-fitting procedures 
would yield a solution to the inverse problem with a guaranteed accuracy. However, 
overwhelming majority of practitioners are using parameter-fitting procedures. In 
dozens of published papers the results obtained by various parameter-fitting pro- 
cedures look quite good. The explanation, in most of the cases is simple: the 
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authors know the answer beforehand, and it is usuaUy not difficult to parametrize 
the unknown function so that the exact solution is well approximated by a function 
from a finite-parametric family, and since the authors know a priori the exact an- 
swer, they may choose numerically the values of the parameters which yield a good 
approximation of the exact solution. When can one rely on the results obtained 
by parameter-fitting procedures? Unfortunately, there is no rigorous and complete 
answer to this question, but some recommendations are given in Section 4- 

In this paper the authors present their recent results which are based on spe- 
cially designed parameter-fitting procedures. Before describing them, let us mention 
that usually in a numerical solution of an inverse scattering problem one uses a reg- 
ularization procedure, e.g. a variational regularization, spectral cut-off, iterative 
regularization, DSM (the dynamical systems method), quasi-solutions, etc, see e.g. 
|78j . j73| . This general theoretical framework is well established in the theory of 
ill-posed problems, of which the inverse scattering problems represent an important 
class. This framework is needed to achieve a stable method for assigning a solu- 
tion to an ill-posed problem, usually set in an infinite dimensional space. The goal 
of this paper is to present optimization algorithms already in a finite dimensional 
setting of a Direct or Inverse scattering problem. 

In Section 2 the problem of finding small subsurface inclusions from surface 
scattering data is investigated f|61p. This (geophysical) Inverse Scattering problem 
is reduced to an optimization problem. This problem is solved by the Hybrid 
Stochastic-Deterministic minimization algorithm ( 27 ). It is based on a genetic 
minimization algorithm ideas for its random (stochastic) part, and a deterministic 
minimization without derivatives used for the local minimization part. 

In Section 3 a similar approach is used to determine layers in a particle sub- 
jected to acoustic or electromagnetic waves. The global minimization algorithm 
uses Rinnooy Kan and Timmer's Multilevel Single-Linkage Method for its stochas- 
tic part. 

In Section 4 we discuss an Inverse potential scattering problem appearing in a 
quantum mechanical description of particle scattering experiments. The central fea- 
ture of the minimization algorithm here is the Stability Index Method ( 72 ). This 
general approach estimates the size of the minimizing sets, and gives a practically 
useful stopping criterion for global minimization algorithms. 

In Section 5 Ramm's method for solving 3D inverse scattering problem with 
fixed-energy data is presented following |79) . see also |K8) . The cases of exact and 
noisy discrete data are considered. Error estimates for the inversion algorithm are 
given in both cases of exact and noisy data. Comparison of the Ramm's inversion 
method with the inversion based on the Dirichlet-to-Neumann map is given and it 
is shown that there are many more numerical difficulties in the latter method than 
in Ramm's method. 

In Section 6 an Obstacle Direct Scattering problem is treated by a novel Mod- 
ified Raylcigh Conjecture (MRC) method. It was introduced in |69| and applied 
in _^30j ^2^ and 7^. MRC's performance is compared favorably to the well known 
Boundary Integral Equation Method, based on the properties of the single and 
double-layer potentials. A special minimization procedure allows us to inexpen- 
sivly compute scattered fields for several 2D and 3D obstacles having smooth as 
well as nonsmooth surfaces. 
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In Section 7 a new Support Function Method (SFM) is used to determine the 
location of an obstacle (cf |48| . |50| . |31p . Unlike other methods, the SFM can 
work with limited data. It can also be used for Inverse scattering problems with 
unknown scattering conditions on its boundary (e.g. soft, or hard scattering). 

Finally, in Section 8, we present an analysis of another popular method for 
Inverse scattering problems, the Linear Sampling Method (LSM), and show that 
both theoretically and computationally the method fails in many aspects. 



2.1. Problem description. In many applications it is desirable to find small 
inhomogeneities from surface scattering data. For example, such a problem arises 
in ultrasound mammography, where small inhomogeneities are cancer cells. Other 
examples include the problem of finding small holes and cracks in metals and other 
materials, or the mine detection. The scattering theory for small scatterers orig- 
inated in the classical works of Lord Rayleigh (1871). Rayleigh understood that 
the basic contribution to the scattered field in the far-field zone comes from the 
dipole radiation, but did not give methods for calculating this radiation. Analyt- 
ical formulas for calculating the polarizability tensors for homogeneous bodies of 
arbitrary shapes were derived in 50 (see also references therein). These formulas 
allow one to calculate the S'-matrix for scattering of acoustic and electromagnetic 
waves by small bodies of arbitrary shapes with arbitrary accuracy. Inverse scatter- 
ing problems for small bodies are considered in |49| and |56| . In |61| the problem 
of identification of small subsurface inhomogeneities from surface data was posed 
and its possible applications were discussed. 

In the context of a geophysical problem, let y e R'^ be a point source of 
monochromatic acoustic waves on the surface of the earth. Let u{x,y,k) be the 
acoustic pressure at a point a: G M'^, and k > he the wavenumber. The governing 
equation for the acoustic wave propagation is: 



where x = {xi, X2, X3), v{x) is the inhomogeneity in the velocity profile, and 
u{x, y, k) satisfies the radiation condition at infinity, i.e. it decays sufficiently fast 
as 00. 

Let us assume that v{x) is a bounded function vanishing outside of the domain 
D = U^^iDm which is the union of AI small nonintersecting domains Dm, all 
of them are located in the lower half-space = {a; : X3 < 0}. Smallness is 
understood in the sense kp <C 1, where p := i maxi<m<M{diamDj„}j and diam 
D is the diameter of the domain D. Practically kp ^ 1 means that kp < 0.1. In 
some cases kp < 0.2 is sufficient for obtaining acceptable numerical results. The 
background velocity in (|2.1|l equals to 1, but we can consider the case of fairly 
general background velocity 56 . 

Denote Zm and Vm the position of the center of gravity of Dm, and the total 
intensity of the m-th inhomogeneity Vm ■= v{x)dx. Assume that Vm 7^ 0. Let 
P be the equation of the surface of the earth: 



2. Identification of small subsurface inclusions. 



(2.1) 



[V^ -I- fc2 -f k^v{x)\ u 



5{x - y) in 



(2.2) 



P — {x = {xi,X2,x:i) e : 



^3 = 0}. 



The inverse problem to be solved is: 
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IP: Given u{x,y, k) for all source- detector pairs {x,y) on P at a fixed k > 0, 
find the number M of small inhomogeneities, the positions Zm of the inhomo- 
geneities, and their intensities Vm- 

Practically, one assumes that a fixed wavenumber fc > 0, and J source-detector 
pairs {xj,yj),j — 1,2,..., J, on P are known together with the acoustic pressure 
measurements ^(xj, j/j, /c). Let 

/r, f i\ exp{ik\x -y\) 

2.3 g{x,y,k) -.^ —— — , x,y e P., 

4:Tr\x — y\ 

(2.4) Gj(z) := G{xj,yj,z) := g{xj,z,k)g{yj,z,k), Xj,yj&P, z e , 



ir,r\ f _ u{xj,yj,k)- g{xj,yj,k) 



and 

,7 

(2.6) ^{zi,...,zm, vi,...,vm) :^ 



M 

fj — ^ Gj{Zm)v„ 
m—1 



The proposed method for solving the (IP) consists of finding the global mini- 
mizer of function H2.5|l . This minimizer (zi, . . . , zm, vi, . . . , u^/) gives the estimates 
of the positions Zm of the small inhomogeneities and their intensities Vm ■ See j61) 
for a justification of this approach. 

The function (f> depends on M unknown points Zm G and M unknown 
parameters Vm, 1 < m < M. The number M of the small inhomogeneities is also 
unknown, and its determination is a part of the minimization problem. 

2.2. Hybrid Stochastic-Deterministic Method. Let the inhomogeneities 
be located within the box 

(2.7) B = {{xi,X2,X3) : —a < xi < a, —b<X2<b, < X3 < c} , 
and their intensities satisfy 

(2.8) 

The box is located above the earth surface for a computational convenience. 

Then, given the location of the points zi, Z2, . . . , zm, the minimum of with 
respect to the intensities vi,V2, ■ ■ ■ , f m can be found by minimizing the resulting 
quadratic function in H2.6() over the region satisfying l|2.8|l . This can be done using 
normal equations for (|2.6|) and projecting the resulting point back onto the region 
defined by l|2.8|l . Denote the result of this minimization by $, that is 

^2 $(zi,Z2, . . . ,zm) = min{$(zi,Z2, . . . ,ZM,fi,W2, • ■ • jI'm) : 

< t)„i < V,nax , 1 < m < M} 

Now the original minimization problem for $(zi, Z2, . . . , zm, vi,V2, ■ ■ ■ ,vm) is 
reduced to the 3 Af -dimensional constrained minimization for <i>(zi, Z2, . . . , zm)'- 

(2.10) $(zi,Z2, . . . ,zm) = min, z^ e B , l<m<M. 

Note, that the dependency of $ on its 3M variables (the coordinates of the 
points Zm) is highly nonlinear. In particular, this dependency is complicated by 
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the computation of the minimum in (|2.9|l and the consequent projection onto the 
admissible set B. Thus, an analytical computation of the gradient of $ is not 
computationally efficient. Accordingly, the Powell's quadratic minimization method 
was used to find local minima. This method uses a special procedure to numerically 
approximate the gradient, and it can be shown to exhibit the same type of quadratic 
convergence as conjugate gradient type methods (see [11) 1. 

In addition, the exact number of the original inhomogeneities Morig is unknown, 
and its estimate is a part of the inverse problem. In the HSD algorithm described 
below this task is accomplished by taking the initial number AI sufficiently large, 
so that 

(2.11) Mo„g<M, 

which, presumably, can be estimated from physical considerations. After all, our 
goal is to find only the strongest inclusions, since the weak ones cannot be distin- 
guished from background noise. The Reduction Procedure (see below) allows the 
algorithm to seek the minimum of $ in a lower dimensional subsets of the admissible 
set B, thus finding the estimated number of inclusions M. Still another difficulty in 
the minimization is a large number of local minima of $. This phenomenon is well 
known for objective functions arising in various inverse problems, and we illustrate 
this point in Figure H 

For example, let Morig = 6, and the coordinates of the inclusions, and their 
intensities (zi, . . . , 5g, -Di, . . . , -Dg) be as in Tabled Figure Q shows the values of the 
function $(zr, Z2, Z3, Z4, Z5, zg), where 

= (r, 0,0.520), -2 <r<2 

and 

Z2 = (-1,0.3,0.580). 
The plot shows multiple local minima and almost flat regions. 

A direct application of a gradient type method to such a function would result 
in finding a local minimum, which may or may not be the sought global one. In 
the example above, such a method would usually be trapped in a local minimum 
located at r = —2, r ~ —1.4, r = —0.6, r = 0.2 or r = 0.9, and the desired 
global minimum at r = 1.6 would be found only for a sufhciently close initial guess 
1.4 < r < 1.9. Various global minimization methods are known (see below), but 
we found that an efficient way to accomplish the minimization task for this Inverse 
Problem was to design a new method (HSD) combining both the stochastic and the 
deterministic approach to the global minimization. Deterministic minimization al- 
gorithms with or without the gradient computation, such as the conjugate gradient 
methods, are known to be efficient (see |11L I20L l46l I38| ). and j82j . However, the 
initial guess should be chosen sufficiently close to the sought minimum. Also such 
algorithms tend to be trapped at a local minimum, which is not necessarily close 
to a global one. A new deterministic method is proposed in [5] and |6], which is 
quite efficient according to [6|. On the other hand, various stochastic minimization 
algorithms, e.g. the simulated annealing method J39' 140' , are more likely to find a 
global minimum, but their convergence can be very slow. We have tried a variety 
of minimization algorithms to find an acceptable minimum of $. Among them 
were the Levenberg-Marquardt Method, Conjugate Gradients, Downhill Simplex, 
and Simulated Annealing Method. None of them produced consistent satisfactory 
results. 
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I 0752^ 1 

-2-1012 

r 

Figure 1. Objective function ^{zr, Z2, S3, Z4, Z5, zg) , ~2 < r < 2 

Among minimization methods combining random and deterministic searches 
we mention Deep's method |19| and a variety of clustering methods |80) . |81| . An 
apphcation of these methods to the particle identification using light scattering is 
described in |84| . The clustering methods are quite robust (that is, they consis- 
tently find global extrema) but, usually, require a significant computational effort. 
One such method is described in the next section on the identification of layers 
in a multilayer particle. The HSD method is a combination of a reduced sample 
random search method with certain ideas from Genetic Algorithms (see e.g. 33 ). 
It is very efficient and seems especially well suited for low dimensional global mini- 
mization. Further research is envisioned to study its properties in more detail, and 
its applicability to other problems. 

The steps of the Hybrid Stochastic-Deterministic (HSD) method are outlined 
below. Let us call a collection of M points ( inclusion's centers) {zi, Z2, zm}, Zi G 
B a configuration Z . Then the minimization problem 1)2.10(1 is the minimization of 
the objective function $ over the set of all configurations. 

For clarity, let Po = 1, Cs = 0.5, €i = 0.25, €d = 0.1, be the same values as the 
ones used in numerical computations in the next section. 



8 



ALEXANDER G. RAMM AND SEMION GUTMAN 



Generate a random configuration Z . Compute the best fit intensities Vi corre- 
sponding to tliis configuration. If Vi > Vmax, then let Vi := Vmax- If Vi < 0, then 
let Vi := 0. If ^(Z) < Pq^s, then this configuration is a preliminary candidate for 
the initial guess of a deterministic minimization method (Step 1). 

Drop the points Zi G Z such that Vi < Vmax^i- That is, the inclusions with 
small intensities are eliminated (Step 2). 

If two points Zk, Zj G Z are too close to each other, then replace them with one 
point of a combined intensity (Step 3). 

After completing steps 2 and 3 we would be left with N < M points zi,Z2,---,zn 
(after a re-indexing) of the original configuration Z. Use this reduced configuration 
Zred as the starting point for the deterministic restraint minimization in the 3A'^ 
dimensional space (Step 4). Let the resulting minimizer be Zred = {zi, zn). If 
the value of the objective function ^(Zred) < then we are done: Zred is the 
sought configuration containing N inclusions. If ^{Zred) > e, then the iterations 
should continue. 

To continue the iteration, randomly generate M — N points in B (Step 5). Add 
them to the reduced configuration Zred- Now we have a new full configuration Z, 
and the iteration process can continue (Step 1). 

This entire iterative process is repeated Umax times, and the best configuration 
is declared to represent the sought inclusions. 

2.3. Hybrid Stochastic-Deterministic (HSD) Method. Let Pq, Tmax, 

nmax, e«) and e be positive numbers. Let a positive integer M be larger 

than the expected number of inclusions. Let N = 0. 

(1) Randomly generate M — N additional points zat+i, . . . , zm € B to obtain 
a full configuration Z = (zi, . . . , zm)- Find the best fit intensities Vi, i = 
1,2, M. If Vi^> Vmax, then let Vi := Vmax- If Vi < 0, then let Vi := 0. 
Compute Ps = ^{zi, Z2 ■ . ■ , zm)- If Ps < Pq^s then go to step 2, otherwise 
repeat step 1. 

(2) Drop all the points with the intensities Vi satisfying Vi < Vmax^i- Now 
only N < M points Zi, Z2 ■ ■ ■ , Zn (re-indexed) remain in the configuration 
Z. 

(3) If any two points Zm , Zn in the above configuration satisfy | — z„ | < edD, 
where D — diam.{B), then eliminate point , change the intensity of point 
Zm to Vm + w„, and assign N :— N — 1. This step is repeated until no 
further reduction in N is possible. Call the resulting reduced configuration 
with points by Zred- 

(4) Run a constrained deterministic minimization of $ in 3A variables, with 
the initial guess Zred- Let the minimizer be Zred = {zi, - - - ,zn)- If 
P = ^{zi, . . . ,zjv) < e, then save this configuration, and go to step 6, 
otherwise let Pq = P, and proceed to the next step 5. 

(5) Keep intact N points z\,..., Zn- If the number of random configurations 
has exceeded Tmax (the maximum number of random tries), then save the 
configuration and go to step 6, otherwise go to step 1, and use these N 
points there. 

(6) Repeat steps 1 through 5 rimax times. 

(7) Find the configuration among the above Umax ones, which gives the small- 
est value to This is the best fit. 



OPTIMIZATION METHODS IN DIRECT AND INVERSE SCATTERING 9 



Table 1. Actual inclusions. 



Inclusions 




X2 


X3 


V 


1 


1.640 


-0.510 


0.520 


1.200 


2 


-1.430 


-0.500 


0.580 


0.500 


3 


1.220 


0.570 


0.370 


0.700 


4 


1.410 


0.230 


0.740 


0.610 


5 


-0.220 


0.470 


0.270 


0.700 


6 


-1.410 


0.230 


0.174 


0.600 



The Powell's minimization method (see for a detailed description) was used 
for the deterministic part, since this method does not need gradient computations, 
and it converges quadratically near quadratically shaped minima. Also, in step 
1, an idea from the Genetic Algorithm's approach 33 is implemented by keeping 
only the strongest representatives of the population, and allowing a mutation for 
the rest. 

2.4. Numerical results. The algorithm was tested on a variety of config- 
urations. Here we present the results of just two typical numerical experiments 
illustrating the performance of the method. In both experiments the box B is 
taken to be 

B = {(xi, X2, xj) : —a < xi < a, —b<X2<b, < < c} , 

with a — 2, 6=1, c — 1. The wavenumber k — 5, and the effective intensities 
are in the range from to 2. The values of the parameters were chosen as follows 

Po = l, T,nax = 1000, = 0.5, e, = 0.25, = 0.1 , e = 10"^ , = 6 

In both cases we searched for the same 6 inhomogeneities with the coordinates 
2^1, 2:2, 2^3 and the intensities v shown in Tabled 

Parameter M was set to 16, thus the only information on the number of in- 
homogeneities given to the algorithm was that their number does not exceed 16. 
This number was chosen to keep the computational time within reasonable limits. 
Still another consideration for the number M is the aim of the algorithm to find 
the presence of the most influential inclusions, rather then all inclusions, which is 
usually impossible in the presence of noise and with the limited amount of data. 

Experiment 1. In this case we used 12 sources and 21 detectors, all on the 
surface X3 = 0. The sources were positioned at {(— 1.667-f 0.667z, — 0.5+1. Oj, 0), i = 
0, 1, . . . , 5, j = 0, 1}, that is 6 each along two lines X2 = —0.5 and X2 = 0.5. The 
detectors were positioned at {{—2 -I- 0.667z, — 1.0 -|- 1.0j,0), i — 0, 1, ... ,6, j = 
0, 1, 2}, that is seven detectors along each of the three hnes X2 = —1, 2:2 = and 
X2 — 1- This corresponds to a mammography search, where the detectors and the 
sources are placed above the search area. The results for noise level i5 = 0.00 are 
shown in Figure [3 and Table [3 The results for noise level i5 = 0.05 are shown in 
Table El 

Experiment 2. In this case we used 8 sources and 22 detectors, all on 
the surface ^3 = 0. The sources were positioned at {(—1.75 + 0.5i, 1.5,0), i = 
0,1,...,7, j = 0,1}, that is all 8 along the line X2 — 1.5. The detectors were 
positioned at {(-2 + OAi, 1.0 + l.Oj, 0), z = 0, 1, . . . , 10, j = 0, 1}, that is eleven 
detectors along each of the two lines X2 = 1 and X2 = 2. This corresponds to a mine 



10 ALEXANDER G. RAMM AND SEMION GUTMAN 

Table 2. Experiment 1. Identified inclusions, no noise, 6 — 0.00. 



Xl 


X2 


X3 


V 


1.640 


-0.510 


0.520 


1.20000 


-1.430 


-0.500 


0.580 


0.50000 


1.220 


0.570 


0.370 


0.70000 


1.410 


0.230 


0.740 


0.61000 


-0.220 


0.470 


0.270 


0.70000 


-1.410 


0.230 


0.174 


0.60000 



Table 3. Experiment 1. Identified inclusions, S = 0.05. 



Xl 


X2 


X3 


V 


1.645 


-0.507 


0.525 


1.24243 


1.215 


0.609 


0.376 


0.67626 


-0.216 


0.465 


0.275 


0.69180 


-1.395 


0.248 


0.177 


0.60747 



Table 4. Experiment 2. Identified inclusions, no noise, S = 0.00. 



Xl 


X2 


X3 


V 


1.656 


-0.409 


0.857 


1.75451 


-1.476 


-0.475 


0.620 


0.48823 


1.209 


0.605 


0.382 


0.60886 


-0.225 


0.469 


0.266 


0.69805 


-1.406 


0.228 


0.159 


0.59372 



Table 5. Experiment 2. Identified inclusions, S ~ 0.05. 



Xl 


X2 


X3 


V 


1.575 


-0.523 


0.735 


1.40827 


-1.628 


-0.447 


0.229 


1.46256 


1.197 


0.785 


0.578 


0.53266 


-0.221 


0.460 


0.231 


0.67803 



search, where the detectors and the sources must be placed outside of the searched 
ground. The results of the identification for noise level S = 0.00 in the data are 
shown in Figure Inland Table 0| The results for noise level S = 0.05 are shown in 
Table m 

In general, the execution times were less than 2 minutes on a 333MIIz PC. 
As it can be seen from the results, the method achieves a perfect identification in 
the Experiment #1 when no noise is present. The identification deteriorates in the 
presence of noise, as well as if the sources and detectors are not located directly 
above the search area. Still the inclusions with the highest intensity and the closest 
ones to the surface are identified, while the deepest and the weakest are lost. This 
can be expected, since their influence on the cost functional is becoming comparable 
with the background noise in the data. 
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Figure 2. Inclusions and Identified objects for subsurface particle 
identification, Experiment 1, ^ = 0.00. 2:3 coordinate is not shown. 



In summary, the proposed method for the identification of small inclusions can 
be used in geophysics, medicine and technology. It can be useful in the development 
of new approaches to ultrasound mammography. It can also be used for localization 
of holes and cracks in metals and other materials, as well as for finding mines from 
surface measurements of acoustic pressure and possibly in other problems of interest 
in various applications. 

The HSD minimization method is a specially designed low-dimensional mini- 
mization method, which is well suited for many inverse type problems. The prob- 
lems do not necessarily have to be within the range of applicability of the Born 
approximation. It is highly desirable to apply HSD method to practical problems 
and to compare its performance with other methods. 
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Figure 3. Inclusions and Identified objects for for subsurface par- 
ticle identification, Experiment 2, S — 0.00. X3 coordinate is not 
shown. 



3. Identification of layers in multilayer particles. 

3.1. Problem Description. Many practical problems require an identifica- 
tion of the internal structure of an object given some measurements on its surface. 
In this section we study such an identification for a multilayered particle illuminated 
by acoustic or electromagnetic plane waves. Thus the problem discussed here is an 
inverse scattering problem. A similar problem for the particle identification from 
the light scattering data is studied in |84) . Our approach is to reduce the inverse 
problem to the best fit to data multidimensional minimization. 

Let D gR^ he the circle of a radius R> 0, 

(3.1) D„, = {xeR'^:r,n-i<\x\<r,„., m = 1,2, . . . , N} 

and Sm — {x £ : \x\ = r,„} for = < ri < ■ ■ ■ < t'jv < R- Suppose 
that a multilayered scatterer in D has a constant refractive index Um in the region 
Dm , m = 1,2, . . . , N . If the scatterer is illuminated by a plane harmonic wave 
then, after the time dependency is eliminated, the total field u{x) — uo{x) + Us{x) 
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satisfies the Helmholtz equation 

(3.2) Au + Zcq^^O, \x\>rN 

where uo{x) = e*'^''^ " is the incident field and a is the unit vector in the direction 
of propagation. The scattered field Us is required to satisfy the radiation condition 
at infinity, see |50j. 

Let ki = k^Tim- We consider the following transmission problem 

(3.3) AUm + k^Ujn =0 X e Dm , 

under the assumption that the fields u,„ and their normal derivatives are continuous 
across the boundaries Sm , m — 1,2, . . . , N . 

In fact, the choice of the boundary conditions on the boundaries Sm depends 
on the physical model under the consideration. The above model may or may not 
be adequate for an electromagnetic or acoustic scattering, since the model may 
require additional parameters (such as the mass density and the compressibility) to 
be accounted for. However, the basic computational approach remains the same. 
For more details on transmission problems, including the questions on the existence 
and the uniqueness of the solutions, see [2], j66, . and ,22, . 

The Inverse Problem to be solved is: 

IPS: Given u{x) for all x G S = {x : \x\ = R) at a fixed fco > 0, find the 
number N of the layers, the location of the layers, and their refractive indices 
Um , m — 1,2, N in (|3.3|l . 

Here the IPS stands for a Single frequency Inverse Problem. Numerical expe- 
rience shows that there are some practical difficulties in the successful resolution of 
the IPS even when no noise is present, see 28 . While there are some results on the 
uniqueness for the IPS (see T, ^6" ) , assuming that the refractive indices are known, 
and only the layers are to be identified, the stability estimates are few, (58j .|59 | . 
j68| . The identification is successful, however, if the scatterer is subjected to a 
probe with plane waves of several frequencies. Thus we state the Multifrequency 
Inverse Problem: 

IPM: Given vP{x) for all x Cz S = {x : \x\ = R) at a finite number P of wave 
numbers k\f^ > 0, find the number N of the layers, the location of the layers, and 
their refractive indices Um , m = 1,2, N in (|3.3() . 

3.2. Best Fit Profiles and Local Minimization Methods. If the refrac- 
tive indices Um are sufficiently close to 1, then we say that the scattering is weak. 
In this case the scattering is described by the Born approximation, and there are 
methods for the solution of the above Inverse Problems. See |15| . |50| and |56| for 
further details. In particular, the Born inversion is an ill-posed problem even if the 
Born approximation is very accurate, see j55| . or |53| . When the assumption of 
the Born approximation is not appropriate, one matches the given observations to 
a set of solutions for the Direct Problem. Since our interest is in the solution of the 
IPS and IPM in the non-Born region of scattering, we choose to follow the best fit 
to data approach. This approach is used widely in a variety of applied problems, 
see e. g. 

Note, that, by the assumption, the scatterer has the rotational symmetry. Thus 
we only need to know the data for one direction of the incident plane wave. For 
this reason we fix a = (1, 0) in H3.2() and define the (complex) functions 

(3.4) g^^P\e), Q<e<2n, p=l,2,...,P, 
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to be the observations measured on the surface S of the ball D for a finite set of 
free space wave numbers kjf^ . 

Fix a positive integer M. Given a configuration 

(3.5) Q = (ri,r2, . . . ,rM,ni,n2, • ■ • ,«A/) 

we solve the Direct Problem (|3.2|) - (|3.3f) (for each free space wave number kjf'^) 
with the layers = {x S : r,m_i < |x| < , m = 1, 2, . . . , M}, and the 
corresponding refractive indices rim, where ro = 0. Let 

(3.6) w<^p\e)^u<^pHx)l^^. 

Fix a set of angles Q — {61,62, ■ ■ ■ ,0l) and let 

L 

(3.7) \\w\\2^{J2w'{6i)r/\ 

1=1 



Define 

(3.8) ^{ri,r2, . ■ . ,rM,ni,n2, . ■ . ,nM) = p X! 



,(P)-g(P)\\l 



where the same set Q is used for g^P^ as for w^P> . 

We solve the IPM by minimizing the above best fit to data functional <& over 
an appropriate set of admissible parameters Aadm C R"^^^ . 

It is reasonable to assume that the underlying physical problem gives some 
estimate for the bounds n/oto and rihigh of the refractive indices rim as well as for 
the bound M of the expected number of layers N. Thus, 

(3.9) 

Aadm C {{ri,r2,. ■ ■ ,rM,ni,n2,. . ■ ,nM) ■ < n < R , niow < n-m < n-Mgh} ■ 
Note, that the admissible configurations must also satisfy 



(3.10) ri < r2 < rg < • • • < TM . 

It is well known that a multidimensional minimization is a difficult problem, 
unless the objective function is "well behaved" . The most important quality of such 
a cooperative function is the presence of just a few local minima. Unfortunately, 
this is, decidedly, not the case in many applied problems, and, in particular, for the 
problem under the consideration. 

To illustrate this point further, let P be the set of three free space wave numbers 
kX chosen to be 



(3.11) P = {3.0, 6.5, 10.0}. 

Figure ^ shows the profile of the functional $ as a function of the variable 
t, 0.1 <t< 0.6 in the configurations qt with 

'0.49 0<|a;|<t 
n{x)^{9.0 t<|2:|<0.6 
1.0 0.6<|2:|<1.0 
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Thus the objective function $ has many local minima even along this arbitrarily 
chosen one dimensional cross-section of the admissible set. There are sharp peaks 
and large gradients. Consequently, the gradient based methods (see ^1, 20, 231 
l34U38ll46] l. would not be successful for a significant portion of this region. It is also 
appropriate to notice that the dependency of $ on its arguments is highly nonlinear. 
Thus, the gradient computations have to be done numerically, which makes them 
computationally expensive. More importantly, the gradient based minimization 
methods (as expected) perform poorly for these problems. 

These complications are avoided by considering conjugate gradient type algo- 
rithms which do not require the knowledge of the derivatives at all, for example the 
Powell's method. Further refinements in the deterministic phase of the minimiza- 
tion algorithm are needed to achieve more consistent performance. They include 
special line minimization, and Reduction procedures similar to the ones discussed 
in a previous section on the identification of underground inclusions. We skip the 
details and refer the reader to j28l . 

In summary, the entire Local Minimization Method (LMM) consists of the 
following: 

Local Minimization Method (LMM). 

(1) Let your starting configuration be Qo — (ri, r2, . . . , tm, ni, n2, . . . , um)- 



16 



ALEXANDER G. RAMM AND SEMION GUTMAN 



(2) Apply the Reduction Procedure to Qo, and obtain a reduced configuration 
Qq containing M'' layers. 

(3) Apply the Basic Minimization Method in Aadm f] '^i^h the starting 
point Qq, and obtain a configuration Qi. 

(4) Apply the Reduction Procedure to Qi, and obtain a final reduced config- 
uration Ql. 



3.3. Global Minimization Methods. Given an initial configuration Qq a lo- 
cal minimization method finds a local minimum near Qo- On the other hand, global 
minimization methods explore the entire admissible set to find a global minimum of 
the objective function. While the local minimization is, usually, deterministic, the 
majority of the global methods are probabilistic in their nature. There is a great 
interest and activity in the development of efficient global minimization methods, 
see e.g. Among them are the simulated annealing method ( 39 , 40 ), var- 

ious genetic algorithms |33j . interval method, TRUST method (0,0), etc. As 
we have already mentioned before, the best fit to data functional $ has many nar- 
row local minima. In this situation it is exceedingly unlikely to get the minima 
points by chance alone. Thus our special interest is for the minimization methods, 
which combine a global search with a local minimization. In j27| we developed 
such a method (the Hybrid Stochastic-Deterministic Method), and applied it for 
the identification of small subsurface particles, provided a set of surface measure- 
ments, see Sections 2.2-2.4. The HSD method could be classified as a variation of 
a genetic algorithm with a local search with reduction. In this paper we consider 
the performance of two algorithms: Deep's Method, and Rinnooy Kan and Tim- 
mer's Multilevel Single-Linkage Method. Both combine a global and a local search 
to determine a global minimum. Recently these methods have been applied to a 
similar problem of the identification of particles from their light scattering charac- 
teristics in |84l . Unlike '84', our experience shows that Deep's method has failed 
consistently for the type of problems we are considering. See |19j and |84| for more 
details on Deep's Method. 

Multilevel Single-Linkage Method (MSLM). Rinnooy Kan and Timmcr in j80| 
and |81j give a detailed description of this algorithm. Zakovic et. al. in |84) 
describe in detail an experience of its application to an inverse light scattering 
problem. They also discuss different stopping criteria for the MSLM. Thus, we 
only give here a shortened and an informal description of this method and of its 
algorithm. 

In a pure Random Search method a batch H of L trial points is generated 
in Aadm using a uniformly distributed random variable. Then a local search is 
started from each of these L points. A local minimum with the smallest value of <& 
is declared to be the global one. 

A refinement of the Random Search is the Reduced Sample Random Search 
method. Here we use only a certain fixed fraction 7 < 1 of the original batch of L 
points to proceed with the local searches. This reduced sample Hred of points 
is chosen to contain the points with the smallest jL values of among the original 
batch. The local searches are started from the points in this reduced sample. 

Since the local searches dominate the computational costs, we would like to 
initiate them only when it is truly necessary. Given a critical distance d we define 
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a cluster to be a group of points located within the distance d of each other. Intu- 
itively, a local search started from the points within a cluster should result in the 
same local minimum, and, therefore, should be initiated only once in each cluster. 

Having tried all the points in the reduced sample we have an information on the 
number of local searches performed and the number of local minima found. This 
information and the critical distance d can be used to determine a statistical level of 
confidence, that all the local minima have been found. The algorithm is terminated 
(a stopping criterion is satisfied) if an a priori level of confidence is reached. 

If, however, the stopping criterion is not satisfied, we perform another iteration 
of the MSLM by generating another batch of L trial points. Then it is combined 
with the previously generated batches to obtain an enlarged batch of jL points 
(at iteration j), which leads to a reduced sample Hj^^^ of ^jL points. According to 
MSLM the critical distance d is reduced to dj, (note that ^ as j — * cxd, since 
we want to find a minimizer), a local minimization is attempted once within each 
cluster, the information on the number of local minimizations performed and the 
local minima found is used to determine if the algorithm should be terminated, etc. 

The following is an adaptation of the MSLM method to the inverse scattering 
problem presented in Section 3.1. The LMM local minimization method introduced 
in the previous Section is used here to perform local searches. 

MSLM. (at iteration j). 

(1) Generate another batch of L trial points (configurations) from a random 
uniform distribution in Aadm- Combine it with the previously generated 
batches to obtain an enlarged batch of jL points. 

(2) Reduce to the reduced sample Hij,^^ of -f jL points, by selecting the 
points with the smallest jjL values of $ in . 

(3) Calculate the critical distance dj by 



2 ; jL 



■) 



rim _^-i/2 ,M\ ji^olnjLy/^ 

« j = TT ' I i I 1 + — 1 [Uhigh - niow) .j^ I 



d 



(4) Order the sample points in Hi^^^ so that ^{Qi) < $(Qi+i)) « = 1, • • • , 'jjL. 
For each value of i, start the local minimization from Qi, unless there exists 
an index k < i, such that \\Qk — Qi\\ < dj. Ascertain if the result is a 
known local minimum. 

(5) Let K be the number of local minimizations performed, and W be the 
number of different local minima found. Let 



K-W-2 
The algorithm is terminated if 

(3.12) Wtot<W + 0.5. 

Here F is the gamma function, and cr is a fixed constant. 
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A related algorithm (the Mode Analysis) is based on a subdivision of the ad- 
missible set into smaller volumes associated with local minima. This algorithm is 
also discussed in '80j and |81) . From the numerical studies presented there, the 
authors deduce their preference for the MSLM. 

The presented MSLM algorithm was successful in the identification of various 
2D layered particles, see |28| for details. 



4. Potential scattering and the Stability Index method. 

4.1. Problem description. Potential scattering problems are important in 
quantum mechanics, where they appear in the context of scattering of particles 
bombarding an atom nucleus. One is interested in reconstructing the scattering 
potential from the results of a scattering experiment. The examples in Section 4 
deal with finding a spherically symmetric {q — q{r), r — \x\) potential from the 
fixed-energy scattering data, which in this case consist of the fixed-energy phase 
shifts. In |60j . ^8 , 78 and j79| the three-dimensional inverse scattering problem 
with fixed-energy data is treated. 

Let q{x), a; G R^, be a real- valued potential with compact support. Let R > 
be a number such that q{x) = for |a::| > R. We also assume that q £ 
L'^{Br) , Br^ {x ■.\x\ < R,x e 9?}. Let 5^ be the unit sphere, and a e S"^. For 
a given energy fc > the scattering solution ip{x, a) is defined as the solution of 

(4.1) AV' + fcV - 9(a;)'0 = , xe«? 
satisfying the following asymptotic condition at infinity: 

(4.2) ^ = ^0 + 1^, ^o:=e*''"-", a & S\ 



(4.3) lim 
It can be shown, that 



|a:|— r 



dv 

— ikv 

or 



ds = 0. 



(4.4) ip(x, a) — + A(a' , a, k) ^o\— ],asr—>oo, — = a' r := \x\. 

r \r J r 

The function A{a' , a, k) is called the scattering amplitude, a and a' are the 
directions of the incident and scattered waves, and k^ is the energy, see !44J , j56| . 

For spherically symmetric scatterers q(x) = q{r) the scattering amplitude sat- 
isfies A{a' , a, fc) = A{a' ■ a,k). The converse is established in |52| . Following |63 |. 
the scattering amplitude for q — q(r) can be written as 



(4.5) A{a',a,k) = J2J2 Ai{k)Yi^{a')Yi,n{a) , 

where Yim are the spherical harmonics, normalized in L^(5^), and the bar denotes 
the complex conjugate. 
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The fixed-energy phase shifts — tt < Si < tt {Si — d{l,k), k > is fixed) are 
related to Ai{k) (see e.g., |63j ) by the formula: 

(4.6) Ai{k) = ^e^'' sm{Si) . 

Several parameter-fitting procedures were proposed for calculating the poten- 
tials from the fixed-energy phase shifts, (by Fiedledey, Lipperheide, Hooshyar and 
Razavy, loannides and Mackintosh, Newton, Sabatier, May and Scheid, Ramm and 
others). These works are referenced and their results are described in |13| and 
Recent works gHl EE EZl ESI and |S3 IMl EH, |M] present new numeri- 
cal methods for solving this problem. In |71| (also see |74l 178) 1 it is proved that 
the R.Newton-P.Sabatier method for solving inverse scattering problem the fixed- 
energy phase shifts as the data (see |44l I13| ) is fundamentally wrong in the sense 
that its foundation is wrong. In j70| a counterexample is given to a uniqueness 
theorem claimed in a modification of the R.Newton's inversion scheme. 

Phase shifts for a spherically symmetric potential can be computed by a variety 
of methods, e.g. by a variable phase method described in |12) . The computation 
involves solving a nonlinear ODE for each phase shift. However, if the potential 
is compactly supported and piecewise-constant, then a much simpler method de- 
scribed in U and |72| can be used. We refer the reader to these papers for details. 



Let qoir) be a spherically symmetric piecewise-constant potential, {S{kJ)} 



N 
1=1 

be the set of its phase shifts for a fixed k > and a sufficiently large N. Let q{r) 
be another potential, and let {6{k,l)}j^^ be the set of its phase shifts. 
The best fit to data function ^{q, k) is defined by 

The phase shifts are known to decay rapidly with see |62| . Thus, for suffi- 
ciently large A^, the function $ is practically the same as the one which would use 
aU the shifts in The inverse problem of the reconstruction of the potential 

from its fixed-energy phase shifts is reduced to the minimization of the objective 
function $ over an appropriate admissible set. 

4.2. Stability Index Minimization Method. Let the minimization prob- 
lem be 

(4.8) min{$(q) : q e Aadm} 

Let go be its global minimizcr. Typically, the structure of the objective function 
$ is quite complicated: this function may have many local minima. Moreover, the 
objective function in a neighborhood of minima can be nearly flat resulting in large 
minimizing sets defined by 



(4.9) S, = {qe Aadm ■■ Hq) < $(9o) + e} 

for an e > 0. 

Given an e > 0, let he the diameter of the minimizing set S^, which we call 
the Stability Index of the minimization problem i4-<^ - 
Its usage is explained below. 



20 



ALEXANDER G. RAMM AND SEMION GUTMAN 



One would expect to obtain stable identification for minimization problems 
with small (relative to the admissible set) stability indices. Minimization problems 
with large stability indices have distinct minimizers with practically the same val- 
ues of the objective function. If no additional information is known, one has an 
uncertainty of the minimizer's choice. The stability index provides a quantitative 
measure of this uncertainty or instability of the minimization. 

If < T], where is an a priori chosen treshold, then one can solve the global 
minimization problem stably. In the above general scheme it is not discussed in 
detail what are possible algorithms for computing the Stability Index. 

One idea to construct such an algorithm is to iteratively estimate stability 
indices of the minimization problem, and, based on this information, to conclude if 
the method has achieved a stable minimum. 

One such algorithm is an Iterative Reduced Random Search (IRRS) method, 
which uses the Stability Index for its stopping criterion. Let a batch H oi L trial 
points be randomly generated in the admissible set Aadm- Let 7 be a certain fixed 
fraction, e.g., 7 = 0.01. Let Smin be the subset of H containing points {pi} with the 
smallest 7L values of the objective function $ in iJ. We call Smin the minimizing 
set. If all the minimizers in Smin are close to each other, then the objective function 
$ is not flat near the global minimum. That is, the method identifies the minimum 
consistently. Let || • || be a norm in the admissible set. 

Let 

e= max - min ^{pj) 

and 

(4.10) = diam{Smin) = max{||pj - pj\\ : pi,pj e Smin} ■ 

Then can be considered an estimate for the Stability Index of the 
minimization problem. The Stability Index reflects the size of the minimizing sets. 
Accordingly, it is used as a self-contained stopping criterion for an iterative min- 
imization procedure. The identification is considered to be stable if the Stability 
Index < t], for an a priori chosen rj > 0. Otherwise, another batch of L trial 
points is generated, and the process is repeated. We used (3 = 1.1 as described be- 
low in the stopping criterion to determine if subsequent iterations do not produce 
a meaningful reduction of the objective function. 

More precisely 

Iterative Reduced Random Search (IRRS). (at the j— th iteration). 
Fix < 7 < 1, /? > 1, 7? > and Nmax- 

(1) Generate another batch of L trial points in Aadm using a random 
distribution. 

(2) Reduce to the reduced sample i^^in 7^ points by selecting the 
points in with the smallest 7L values of <&. 

(3) Combine H^^^^ with obtained at the previous iteration. Let S'^„j„ 
be the set of 7L points from Hi,- U Hi'~- with the smallest values of 
(Use for i = 1). 

(4) Compute the Stability Index (diameter) of by = max{||pj — 

Pfell : Pi,Pk G Smin} ■ 

(5) Stopping criterion. 
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Let p e S'^i„ be the point with the smallest value of $ in S^^^ (the 
global minimizer). 

If < rj, then stop. The global minimizer is p. The minimization is 
stable. 

If > T] and ^{q) < (3^{p) ■ q G S^^^^, then stop. The minimization 
is unstable. The Stability Index is the measure of the instability of 
the minimization. 

Otherwise, return to step 1 and do another iteration, unless the max- 
imum number of iterations Nmax is exceeded. 

One can make the stopping criterion more meaningful by computing a normal- 
ized stability index. This can be achieved by dividing by a fixed normalization 
constant, such as the diameter of the entire admissible set Aadm- To improve the 
performance of the algorithm in specific problems we found it useful to modify 
(IRRS) by combining the stochastic (global) search with a deterministic local min- 
imization. Such Hybrid Stochastic-Deterministic (HSD) approach has proved to be 
successful for a variety of problems in inverse quantum scattering (see |72L I28LI67] ) 
as well as in other applications (see I27LI26| ). A somewhat different implementation 
of the Stability Index Method is described in 2Q . 

We seek the potentials q{r) in the class of piecewise-constant, spherically sym- 
metric real-valued functions. Let the admissible set be 
(4.11) 

Aadm C {(ri, r2, . . . , TM, gi, 92, ■ ■ • , qhl) ■ <ri < R, qiow < q-m < qhigh} , 

where the bounds qiow and qhigh for the potentials, as well as the bound M on 
the expected number of layers are assumed to be known. 

A configuration (ri, r2, . . . , ri\i, 51, (72, • • • , 9a/) corresponds to the potential 

(4.12) q{r) = , for r™_i < r < r„ , 1 < m < Af , 

where tq — and q{r) = for r > tm — R. 

Note, that the admissible configurations must also satisfy 

(4.13) ri < r2 < < • • • < TM . 

We used /3 = 1.1, e = 0.02 and jmax ~ 30. The choice of these and other 
parameters (L = 5000, 7 — 0.01, v — 0.16 ) is dictated by their meaning in the 
algorithm and the comparative performance of the program at their different values. 
As usual, some adjustment of the parameters, stopping criteria, etc., is needed to 
achieve the optimal performance of the algorithm. The deterministic part of the 
IRRs algorithm was based on the Powell's minimization method, one-dimensional 
minimization, and a Reduction procedure similar to ones described in the previous 
section 3, see |72| for details. 

4.3. Numerical Results. We studied the performance of the algorithm for 
3 different potentials qi{r), z = 1, 2, 3 chosen from the physical considerations. 

The potential qsir) — —10 for < r < 8.0 and 93 = for r > 8.0 and a wave 
number fc = 1 constitute a typical example for elastic scattering of neutral particles 
in nuclear and atomic physics. In nuclear physics one measures the length in units 
of fm — lO^^^m, the quantity 93 in units of 1/fm^, and the wave number in units of 
1/fm. The physical potential and incident energy are given by V{r) — ^qsir) and 
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E = respectively, here h := ^ , h = 6.62510 erg-s is the Planck constant, 

he — 197.32 MeV-fm, c — 3 ■ 10^ m/sec is the velocity of light, and /i is the mass 
of a neutron. By choosing the mass /x to be equal to the mass of a neutron /i = 
939.6 MeV/c^, the potential and energy have the values of V{r) = -207.2 MeV for 
< ?' < 8.0 fm and E{k =l/fm ) = 20.72 MeV. In atomic physics one uses atomic 
units with the Bohr radius Oq = 0.529 • 10~^°m as the unit of length. Here, r, k 
and 93 are measured in units of ao,l/ao and I/oq, respectively. By assuming a 
scattering of an electron with mass mo — 0.511 MeV/c^, we obtain the potential 
and energy as follows: V{r) = —136 eV for < r < 8ao = 4.23 • 10~^°m and 
E(k = 1/ao) = 13.6 eV. These numbers give motivation for the choice of examples 
applicable in nuclear and atomic physics. 

The method used here deals with finite-range (compactly supported) potentials. 
One can use this method for potentials with the Coulomb tail or other potentials 
of interest in physics, which are not of finite range. This is done by using the 
phase shifts transformation method which allows one to transform the phase shifts 
corresponding to a potential, not of finite range, whose behavior is known for r > a, 
where a is some radius, into the phase shifts corresponding to a potential of finite 
range a (see [2], p. 156). 

In practice differential cross section is measured at various angles, and from 
it the fixed-energy phase shifts are calculated by a parameter-fitting procedure. 
Therefore, we plan in the future work to generalize the stability index method to 
the case when the original data are the values of the differential cross section, rather 
than the phase shifts. 

For the physical reasons discussed above, we choose the following three poten- 
tials: 



In each case the following values of the parameters have been used. The radius 
R of the support of each qi was chosen to be i? = 10.0. The admissible set Aadjn 
(I4.11|l was defined with M = 2. The Reduced Random Search parameters: L = 
5000, 7 0.01, ly = 0.16, e = 0.02, /3 = I.IO, j^ax = 30. The value = 0.1 
was used in the Reduction Procedure during the local minimization phase. The 
initial configurations were generated using a random number generator with seeds 
determined by the system time. A typical run time was about 10 minutes on a 333 
MHz PC, depending on the number of iterations in IRRS. The number N of the 
shifts used in (|4.7|) for the formation of the objective function ^(g) was 31 for all 
the wave numbers. It can be seen that the shifts for the potential decay rapidly 
for fc = 1, but they remain large for k = A. The upper and lower bounds for the 
potentials qiow = —20.0 and qhigh = 0.0 used in the definition of the admissible set 
Aadm were chosen to reflect a priori information about the potentials. 






10.0 



< r < 8.0 



r > 8.0 
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Table 6. Stability Indices for qi{r) identification at different noise 
levels h. 



k 


Iteration 


h = 0.00 


h = 0.01 


h = 0.10 


1.00 


1 


1.256985 


0.592597 


1.953778 




2 


0.538440 


0.133685 


0.799142 




3 


0.538253 


0.007360 


0.596742 




4 


0.014616 




0.123247 




5 






0.015899 


2.00 


1 


0.000000 


0.020204 


0.009607 


2.50 


1 


0.000000 


0.014553 


0.046275 


3.00 


1 


0.000000 


0.000501 


0.096444 


4.00 


1 


0.000000 


0.022935 


0.027214 



The identification was attempted with 3 different noise levels h. The levels are 
h = 0.00 (no noise), h = 0.01 and h — 0.1. More precisely, the noisy phase shifts 
6h{k, I) were obtained from the exact phase shifts S{k, I) by the formula 



6h{k,l) ^ S{k,l){l + {0.5 - z) ■ h) , 

where z is the uniformly distributed on [0, 1] random variable. 

The distance d{pi{r),p2{r)) for potentials in step 5 of the IRRS algorithm was 
computed as 



d{pi{r),p2{r)) = \\pi{r) - P2{r)\\ 

where the norm is the L2-norm in R^. 

The results of the identification algorithm (the Stability Indices) for different 
iterations of the IRRS algorithm are shown in Tables EE 

For example. Table |H1 shows that for k = 2.5, h = 0.00 the Stability Index has 
reached the value 0.013621 after 2 iteration. According to the Stopping criterion for 
IRRS, the program has been terminated with the conclusion that the identification 
was stable. In this case the potential identified by the program was 



-10.000024 < r < 7.999994 
0.0 r > 7.999994 



p{r) 

which is very close to the original potential 



-10.0 < r < 8.0 
0.0 r>8.0 



On the other hand, when the phase shifts of (73 (r) were corrupted by a 10% 
noise (k = 2.5, h = 0.10), the program was terminated (according to the Stopping 
criterion) after 4 iterations with the Stability Index at 0.079241. Since the Stability 
Index is greater than the a priori chosen threshold of e = 0.02 the conclusion is 
that the identification is unstable. A closer look into this situation reveals that the 
values of the objective function ^(pi), pi G Smin (there are 8 elements in Smin) 
are between 0.0992806 and 0.100320. Since we chose (3 = 1.1 the values are within 
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Table 7. Stability Indices for q2{r) identification at different 
levels h. 



k 


TtPTfl.tiOTI. 


h = 0.00 


h = 0.01 


h = 0.10 


1.00 


1 


n 774'^7fi 


598471 


1 08902 




2 


0.773718 


1.027345 


0.023206 




3 


026492 


025593 


023206 




4 


020522 


029533 


024081 




5 


020524 


029533 


024081 




5 


n 000745 


029533 






7 

8 

O 




029533 
099533 






9 




0.029533 






10 




029533 






11 




029619 






12 




025816 






13 




025816 






14 




0.008901 




2.00 


1 


863796 


799356 


981239 




2 


861842 


799356 


029445 




3 


008653 


000993 


029445 




4 






029445 




5 






026513 




6 






0.026513 




7 






0.024881 


2.50 


1 


\ 848910 


1 632298 


894087 




2 


1.197131 


1 632298 


507953 




3 


580361 


1 183455 

-L . -L OOTTtJ O 


025454 




4 


n 03051 6 


528979 






5 


016195 


032661 




3 00 


1 


1 844702 


1 849016 


1 708201 

J- • 1 UO^U J- 




2 


1 649700 


1 789775 


1 519891 




Q 
o 


J. .'4;(JUVJZiU 


A.. 1 OZi 1 1 o 






4 


1 41 02^*^ 


1 4^7020 


1 1 56964 






n 694358 


961 963 


1 1 56964 




U 


n 699080 


961 963 


909681 




7 


u.oyzuoU 


U.yDiZDO 


u.yuzooi 




8 


0.345804 


0.291611 


0.902474 




9 


0.345804 


0.286390 


0.159221 




10 


0.345804 


0.260693 


0.154829 




11 


0.043845 


0.260693 


0.154829 




12 


0.043845 


0.260693 


0.135537 




13 


0.043845 


0.260693 


0.135537 




14 


0.043845 


0.260693 


0.135537 




15 


0.042080 


0.157024 


0.107548 




16 


0.042080 


0.157024 






17 


0.042080 


0.157024 






18 


0.000429 


0.157024 






19 




0.022988 





4.00 



1 0.000000 0.000674 0.050705 
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Table 8. Stability Indices for qair) identification at different noise 
levels h. 



k 


ItcTCbtion 


h = 0.00 


h = 0.01 


h = 0.10 


1.00 


1 


5641 68 


594314 


764340 




2 


024441 


028558 


081888 




3 


024441 


014468 


050755 

\j t\jyj\j 1 yjyj 




4 


024684 








5 


0.024684 








(j 


O.OO-kSOO 






2.00 


1 


684053 


1 450148 


485783 




2 


423283 


792431 


078716 




3 


0.006291 


0.457650 


0.078716 




4 




023157 


078716 




5 






078716 




6 






0.078716 




7 






078716 

\j * yj t yj t i Kj 




8 






078716 




9 






0.078716 




10 






0.078716 




11 






0.078716 


2.50 


1 


0.126528 


993192 


996519 




2 


013621 


105537 


855049 




3 




033694 


849123 




4 




0.020811 


0.079241 


3.00 


1 


0.962483 


1.541714 


0.731315 




2 


0.222880 


0.164744 


0.731315 




3 


0.158809 


0.021775 


0.072009 




4 


0.021366 








5 


0.021366 








() 


().()()141() 






4.00 


1 


1.714951 


1.413549 


0.788434 




2 


0.033024 


0.075503 


0.024482 




3 


0.018250 


0.029385 






4 




0.029421 






5 




0.029421 






() 




U. 015940 





the required 10% of each other. The actual potentials for which the normaUzed 
distance is equal to the Stability Index 0.079241 are 

' -9.997164 < r < 7.932678 
Pi{r) = { -7.487082 7.932678 <r< 8.025500 
0.0 r > 8.025500 

-9.999565 < r < 7.987208 
P2(r) = I -1.236253 7.987208 <r< 8.102628 
0.0 r > 8.102628 



and 
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Table 9. Stability Indices for q2{r) identification for different vaf 
ues of M. 



Iteration 



M =1 M = 2 M = 3 



M = 4 



1 0.472661 1.068993 1.139720 

2 0.000000 0.400304 0.733490 

3 0.000426 0.125855 

4 0.125855 

5 0.033173 

6 0.033173 

7 0.033123 

8 0.000324 
9 



f 
f 












,453076 
,453076 
,899401 
,846117 
,941282 
,655669 
,655669 
,948816 
,025433 
025433 
012586 



10 
11 



with $(pi) = 0.0992806 and $(^2) = 0.0997561. One may conclude from this 
example that the threshold e = 0.02 is too tight and can be relaxed, if the above 
uncertainty is acceptable. 

Finally, we studied the dependency of the Stability Index from the dimension 
of the admissible set Aadm, see (|4.11() . This dimension is equal to 2M , where 
M is the assumed number of layers in the potential. More precisely, M = 3, for 
example, means that the search is conducted in the class of potentials having 3 or 
less layers. The experiments were conducted for the identification of the original 
potential q2(r) with k = 2.0 and no noise present in the data. The results are 
shown in Table |51 Since the potential q2 consists of only one layer, the smallest 
Stability Indices are obtained for M = 1. They gradually increase with M. Note, 
that the algorithm conducts the global search using random variables, so the actual 
values of the indices are different in every run. Still the results show the successful 
identification (in this case) for the entire range of the a priori chosen parameter 
M. This agrees with the theoretical consideration according to which the Stability 
Index corresponding to an ill-posed problem in an infinite-dimensional space should 
be large. Reducing the original ill-posed problem to a one in a space of much lower 
dimension regularizes the original problem. 



5.1. Problem description. In this Section we continue a discussion of the 
Inverse potential scattering with a presentation of Ramm's method for solving in- 
verse scattering problem with fixed-energy data, see |79| . The method is applicable 
to both exact and noisy data. Error estimates for this method are also given. An 
inversion method using the Dirichlet-to-Neumann (DN) map is discussed, the dif- 
ficulties of its numerical implementation are pointed out and compared with the 
difficulties of the implementation of the Ramm's inversion method. See the previous 
Section on the potential scattering for the problem set up. 

5.2. Ramm's inversion method for exact data. The results we describe 
in this Section are taken from (56) and [68 . Assume q Q Q := QaO L°° {M.'^), where 



Qa := {q ■■ q{x) = qix), q{x) G L^{Ba), q{x) = if \x\ > a}, Ba {x : \x\ < 



5. Inverse scattering problem with fixed-energy data. 
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a}. Let A(a',a) be the corresponding scattering amplitude at a fixed energy /c^, 
A; = 1 is taken without loss of generality. One has: 

oo ^ 

(5.1) A{a',a)^Y Ae{a)Ye{a'), Af(a) / A{a' , a)Yi(^da' , 

where is the unit sphere in M"^, Yi{a') = Yi^„i{a'),—t < m < £, are the nor- 
malized spherical harmonics, summation over m is understood in H5.1|) and in (|5.8|) 
below. Define the following algebraic variety: 

3 

(5.2) M := {0 -.0 eC^,e -0 ^1}, e-w:^Y0jWj. 

This variety is non-compact, intersects M.^ over S'^, and, given any ^ e K'^, there 
exist (many) 9, 0' G M such that 

(5.3) 0'-9 = ^, |6'|^oo, 9,9' eM. 

In particular, if one chooses the coordinate system in which ^ = tea, i > 0, 63 is 
the unit vector along the xa-axis, then the vectors 

(5.4) 0' = ^e3 + C2e2 + Ciei, ^ = -^63 + C2e2 + Ciei, Ci+Cl = l-J, 

satisfy H5.3|) for any complex numbers and (2 satisfying the last equation (|5.4|) 
and such that |Ci P + IC2 P 00. There are infinitely many such Ci , C2 € C. Consider 
a subset M' C M consisting of the vectors = (sin?? cos (p, sini? sin cos ■(?), where 
?? and (p run through the whole complex plane. Clearly 9 € M, but M' is a proper 
subset of M. Indeed, any 9 € M with 6*3 ±1 is an element of M' . If ^3 = ±1, 
then cos I? — ±1, so s'm-d = and one gets = (0, 0, ±1) G M'. However, there are 
vectors 9 = {0i,02,l) S M which do not belong to M' . Such vectors one obtains 
choosing 6*1,^2 G C such that 9\ + 02 = 0. There are infinitely many such vectors. 
The same is true for vectors (6*1, 6'2, — 1). Note that in (|5.3|l one can replace M by 
M' for any ^ e RS, ^ ^ 2ez. 

Let us state two estimates proved in j56) : 



(5.5) max|A.(.)|<c(^)^g; 
where c > is a constant depending on the norm HfzllLS^^^), and 



1 r\Ime\ 

(5.6) \Y,{0)\ < ^ I . . , Vr>0, 0&M', 

V47r |j£(r-)| 



where 

(5.7) Mr) := (J;) ' J,+ . M = (i)' [1 + «(!)] ^ ^ 

and J£(r) is the Bessel function regular at r = 0. Note that y^(a'), defined above, 
admits a natural analytic continuation from 5*^ to M by taking -d and (p to be 
arbitrary complex numbers. The resulting 0' S M' C M. 

The series (|5.1f) converges absolutely and uniformly on the sets S*^ x Mc, where 
Mc is any compact subset of M. 
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Fix any numbers ai and b, such that a < ai < b. Let || • || denote the L^(ai < 
|a;| < 6)-norm. If |a;| > a, then the scattering solution is given analytically: 

oo 

(5.8) u{x,a) ^ e"^-"" + J2Ma)Ma')heir), r := \x\ > a, 

where Ai{a) and Yi{a') are defined above, 

h,ir) :=^e'*(^+^/f<,(r), 
V 2r *+2 

H^\r) is the Hankel function, and the normalizing factor is chosen so that hi{r) — 
^[1 + o(l)] as r ^ cx). Define 

(5.9) p{x) -.^ p{x;i^) -.^ e-''^-^ [ u{x,ay{a,9)da - 1, i' e L^{S^). 

Consider the minimization problem 

(5.10) IIpII = inf := d(6l), 

where the infimum is taken over all i/ e i^(5^), and (|5.3(l holds. 
It is proved in |56| that 

(5.11) d(6l) < c|6'|-Mf 61 e M, |6I|>1. 

The symbol \9\ ^ 1 means that \9\ is sufficiently large. The constant c > in (|5.11() 
depends on the norm ||q||L2(B„) but not on the potential q{x) itself. 

An algorithm for computing a function z/(a, 9), which can be used for inversion 
of the exact, fixed-energy, three-dimensional scattering data, is as follows: 

a) Find an approximate solution to (|5.10(l in the sense 

(5.12) \\pix,iy)\\<2di9), 

where in place of the factor 2 in (|5.12() one could put any fixed constant greater 
than 1. 

b) Any such i^{a, 9) generates an estimate of g(^) with the error O ■, \0\ 
oo. This estimate is calculated by the formula 

(5.13) q:=-47r/ A{9' ,a)i^{a,9)da, 

where v{a,9) £ L'^{S'^) is any function satisfying H5.12|l . 
Our basic result is: 



Theorem 5.1. Let ((OJ and (fsT^ hold. Th. 



en 



c 



(5.14) snp\q-qm<T7rr 1^1 

The constant c > in (|5.14|) depends on a norm of q, hut not on a particular q. 

The norm of q in the above Theorem can be any norm such that the set {q : 
\\q\\ < const} is a compact set in L°°[Ba). 

In |56| and j68| an inversion algorithm is formulated also for noisy data, and 
the error estimate for this algorithm is obtained. Let us describe these results. 
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Assume that the scattering data are given with some error: a function As{a' , a) 
is given such that 

(5.15) sup \A{a',a)-As{a',a)\<S. 

We emphasize that As {a', a) is not necessarily a scattering amplitude corre- 
sponding to some potential, it is an arbitrary function in L°°{S'^ x 5^) satisfying 
H5.15|l . It is assumed that the unknown function A(a', a) is the scattering amplitude 
corresponding to a q G Q. 

The problem is: Find an algorithm for calculating qs such that 

(5.16) sup \qs-m\<v{3), 7?(<5) ^ as 5 -> 0, 



and estimate the rate at which ri{5) tends to zero. 

An algorithm for inversion of noisy data will now be described. 
Let 

(5.17) N{S):=^ 



where [a;] is the integer nearest to x > 0, 

NiS) 



ln|ln(5|_ 

(5.18) As{0',a):=y Asi{a)Yeie'), Asi{a) -.^ f Ag^a' , a)Y^da' , 

e=o -Is- 

N{S) 

(5.19) us{x, a) := e*"'" + ^ Asiia)Ye{a')he{r), 

(5.20) ps{x] v) := e"'^-^ / us{x, a)v{a)da - 1, 61 e M, 

(5.21) ^l[5) -.^ e-'^^^^\ 7 = ln— >0, 

a 

(5.22) a{v):=ML2is-), k := |/m0|. 
Consider the variational problem with constraints: 

(5.23) \e\ = sup := d{5), 

(5.24) \B\[\\ps{v)\\+a{v)e''^^i{5)\<c, 9 € M, \0\ = snp ~ ^(S), 

the norm is defined above 1)5. 8|l . and it is assumed that 1)5.3(1 holds, where ^ e M.^ is 
an arbitrary fixed vector, c > is a sufficiently large constant, and the supremum 
is taken over 9 G M and v £ L'^{S'^) under the constraint 15.24|l . By c we denote 
various positive constants. 

Given ^ g one can always find 9 and 9' such that H5.3I) holds. We prove that 
i9((5) — !■ oo, more precisely: 

(5.25) '^(^) ^ ^(hTfli^' 
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Let the pair 9{S) and h's{a,9) be any approximate solution to problem H5.23(l - 
(|5.24(l in the sense that 

(5.26) m\ > 

Calculate 



(5.27) qs:=-4TT Asi9' ,a),^s{a,e)d: 



la. 



Theorem 5.2. // (|5.3|l and (|5.26|l hold, then 
(5.28) sup |$5-q(e)| <cMi^^ as,5-.0, 

where c > is a constant depending on a norm of q. 

In 56 estimates (|5.14|l and (|5.28|l were formulated with the supremum taken 
over an arbitrary large but fixed ball of radius ^o- Here these estimates are im- 
proved: ^0 ^ CO- The key point is: the constant c > in the estimate (jS.llll does 
not depend on 9. 

Remark. In |6Uj (see also |54| and |68p an analysis of the approach to ISP, 
based on the recovery of the DN (Dirichle-to-Neumann) map from the fixed-energy 
scattering data, is given. This approach is discussed below. 

The basic numerical difficulty of the approach described in Theorems 15.11 and 
15.21 comes from solving problems H5.10|l for exact data, and problem (I5.23|l - I|5.24(l 
for noisy data. Solving H5.10|l amounts to finding a global minimizer of a quadratic 
form of the variables ce, if one takes ly in H5.9|l as a linear combination of the 
spherical harmonics: — '^(^q ceYe{a) . If one uses the necessary condition for 
a minimizer of a quadratic form, that is, a linear system, then the matrix of this 
system is ill-conditioned for large L. This causes the main difficulty in the numerical 
solution of (|5.10() . On the other hand, there are methods for global minimization 
of the quadratic functionals, based on the gradient descent, which may be more 
efficient than using the above necessary condition. 

5.3. Discussion of the inversion method which uses the DN map. In 

|60 | the following inversion method is discussed: 



(5.29) q{0= lim / exp(-i6'' • s)(A - Ao)^ds, 

where (|5.3|l is assumed, A is the Dirichlet-to-Neumann (DN) map, t/j is found from 
the equation: 

ie-s 



(5.30) ij{s)=TPo{s)- Gis~t)Bijdt, B:=A-Ao, ^o(s) := e' 

Js 

and G is defined by the formula: 

/r- q1^ r<f \ fa \ ^ f exp(iC • x)d£_ 

(5.31) G(.) = exp(z0 . .)^ J^^ ^TW^- 

The DN map is constructed from the fixed-energy scattering data A{a' , a) by the 
method of [50] (see also 

Namely, given A{a' , a) for all a', a G S^, one finds A using the following steps. 



OPTIMIZATION METHODS IN DIRECT AND INVERSE SCATTERING 31 

Let / e H^/^{S) be given, S* is a sphere of radius a centered at the origin, fi 
are its Fourier coefficients in the basis of the spherical harmonics, 

(5.32) n^^f^f.Y.ix^)!^, r>a, x° -, r := \x\. 
Let 

(5.33) w= g{x, s)a{s)ds, 



where a is some function, which we find below, and g is the Green function (re- 
solvent kernel) of the Schroedinger operator, satisfying the radiation condition at 
infinity. Then 

(5.34) u;+ = + a, 

where N is the outer normal to S, so N is directed along the radius- vector. We 
require w = f on S. Then w is given by 1)5.32(1 in the exterior of S, and 

(5.35) ^-^^Y.f,Y,ix")^ 

1=0 > 

By formulas H5.34|) and H5.35|l . finding A is equivalent to finding a. By H5.33|l . 
asymptotics of w as r := \x\ oo, x/\x\ :— a;*^, is (cf |56| . p. 67): 

(5.36) «; = ^^%^+o(i), 

r Att r 

where u is the scattering solution, 

oo 

(5.37) u(y, -x°) = e-"°-^ + J] M-x'')Y,{y'')h,{\y\). 

^From H5.32|l . H5.36|) and H5.37|l one gets an equation for finding a ('6U', eq. (23), 
see also [56], p. 199): 

(5.38) 1^)^^ ldsais)iuis,-(3),Yim^,^s,^, 
which can be written as a hnear system: 

for the Fourier coefficients of a. The coefficients 

Ai'i := {{A{a' ,a),Yi{a'))L2(s2-^,Yi{a))L2(s2) 

are the Fourier coefficients of the scattering amplitude. Problems (I5.38|) and (|5.39|) 
are very ill-posed (see |6U| for details). 
This approach faces many difficulties: 

1) The construction of the DN map from the scattering data is a very ill-posed 
problem, 

2) The construction of the potential from the DN map is a very difficult problem 
numerically, because one has to solve a Fredholm-type integral equation ( equation 
((5.30|l ) whose kernel contains G, defined in ((5.31|l . This G is a tempered distribu- 
tion, and it is very difficult to compute it. 
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3) One has to calculate a limit of an integral whose integrand grows exponen- 
tially to infinity if a factor in the integrand is not known exactly. The solution 
of equation H5.30|l is one of the factors in the integrand. It cannot be known ex- 
actly in practice because it cannot be calculated with arbitrary accuracy even if the 
scattering data are known exactly. Therefore the limit in formula H5.29|l cannot be 
calculated accurately. 

No error estimates are obtained for this approach. 

In contrast, in Ramm's method, there is no need to compute G, to solve equa- 
tion H5.30|l . to calculate the DN map from the scattering data, and to compute 
the limit H5.29|l . The basic difficulty in Ramm's inversion method for exact data 
is to minimize the quadratic form (|5.1()() . and for noisy data to solve optimization 
problem 1)5. 23|I - H5. 24(1 . The error estimates are obtained for the Ramm's method. 

6. Obstacle scattering by the Modified Rayleigh Conjecture (MRC) 

method. 

6.1. Problem description. In this section we present a novel numerical 
method for Direct Obstacle Scattering Problems based on the Modified Rayleigh 
Conjecture (MRC). The basic theoretical foundation of the method was developed 
in 69 . The MRC has the appeal of an easy implementation for obstacles of com- 
plicated geometry, e.g. having edges and corners. A special version of the MRC 
method was used in 32 to compute the scattered field for 3D obstacles. In our nu- 
merical experiments the method has shown itself to be a competitive alternative to 
the BIEM (boundary integral equations method), see j30| . Also, unlike the BIEM, 
one can apply the algorithm to different obstacles with very little additional effort. 

We formulate the obstacle scattering problem in a 3D setting with the Dirichlet 
boundary condition, but the discussed method can also be used for the Neumann 
and Robin boundary conditions. 

Consider a bounded domain D C M.^, with a boundary S which is assumed to 
be Lipschitz continuous. Denote the exterior domain by D' = M.^\D. Let a, a' G 5^ 
be unit vectors, and be the unit sphere in R^. 

The acoustic wave scattering problem by a soft obstacle D consists in finding 
the (unique) solution to the problem H6.1|l - I(6.2|l : 

(6.1) (V^ + k^)u = in £>', u = on S, 



(6.2) u = uq -t- A{a' , a) ho-, r := \x\ ^ co, a' := — . 

r \r J r 

Here uq := e*'^"'^ is the incident field, v :— u — uq is, the scattered field, A{a' ,a) 
is called the scattering amplitude, its k-dependence is not shown. A: > is the 
wavenumber. Denote 

(6.3) Ae{a):^ [ A{a' ,a)Ye{a')da' , 

where Ye(a) are the orthonormal spherical harmonics. Ye — Yem,—( < m < £. 
Let hi{r) be the spherical Hankcl functions, normalized so that he{r) ~ - — as 
r +00. 
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Informally, the Random Multi-point MRC algorithm can be described as fol- 
lows. 

Fix a J > 0. Let Xj,j = 1, 2, J be a batch of points randomly chosen inside 
the obstacle D. For x & D' , let 

(6.4) "' = r5^' Mx,xj) =Yi{a')he{k\x-Xj\). 

\X Xj\ 

Let g{x) = uo{x), x & S, and minimize the discrepancy 

J L 

(6.5) $(c) = \\g{x) + Y,Y.Ci,jMx,Xj)\\LHs), 

i=i e=o 

over c G , where c = {cej}. That is, the total field u = g{x) -I- w is desired to 
be as close to zero as possible at the boundary S, to satisfy the required condition 
for the soft scattering. If the resulting residual r™™ = min^ is smaller than the 
prescribed tolerance e, than the procedure is finished, and the sought scattered field 
is 

J L 

V^{x) ^'^'^Cl^jijfXx^Xj), XGD', 
3 = 1 £=0 

(see Lemma [6.11 below') . 

If, on the other hand, the residual r™™ > e, then we continue by trying to 
improve on the already obtained fit in l|6.5|l . Adjust the field on the boundary by 
letting g{x) := g{x) + v^{x), x € S. Create another batch of J points randomly 
chosen in the interior of D, and minimize (|6.5|1 with this new g{x). Continue with 
the iterations until the required tolerance e on the boundary S is attained, at the 
same time keeping the track of the changing field Ve- 

Note, that the minimization in H6.5|l is always done over the same number of 
points J. However, the points Xj are sought to be different in each iteration to 
assure that the minimal values of $ are decreasing in consequent iterations. Thus, 
computationally, the size of the minimization problem remains the same. This is 
the new feature of the Random multi-point MRC method, which allows it to solve 
scattering problems untreatable by previously developed MRC methods, see |30| . 

Here is the precise description of the algorithm. 

Random Multi-point MRC. 

For Xj E D, and £>0 functions tpg{x,Xj) are defined as in (|6.4|l . 

(1) Initialization. Fix e > 0, L > 0, J > 0, Nmax > 0. Let n = 0, = 
and g{x) = uq{x), x £ S. 

(2) Iteration. 

(a) Let n := n + 1. Randomly choose J points Xj d D, j — 1,2, . . . , J. 

(b) Minimize 

,7 L 

*(c) = ||5(a:) + ^^C£jV£(2^,2:j)||L2(s) 

over c £ , where c = {cej}- 

Let the minimal value of $ be r"*™. 
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(c) Let 

J L 

Ve{x) := Ve{x) + ^ ^ Cgjlpeix, Xj), X G D'. 
3 = 1 i=0 

(3) Stopping criterion. 

(a) If r""' < e, then stop. 

(b) Ifr™">e, and n ^ iV„a„ let 

J L 

5(2;) := ff(a;) + ^^C£jVf(2::2:j), a; € 5' 
j=i e=Q 

and repeat the iterative step (jSJ. 

(c) If r™™ > e, and n = Nmax, then the procedure failed. 

6.2. Direct scattering problems and the Rayleigh conjecture. Let a 

ball Bfi := {a; : \x\ < R} contain the obstacle D. In the region r > R the solution 
to inU-iniS) is: 

(6.6) u(a;,a) = e''="-"+^A,(a)V^,, ^Pe := Yi{a')he{kr), r > R, a' = ^, 

i=0 

where the sum includes the summation with respect to m, —£ < m < and Aila) 
are defined in 1)6. 

The Rayleigh conjecture (RC) is: the series i)6'.6|) converges up to the boundary 
S (originally RC dealt with periodic structures, gratings). This conjecture is false 
for many obstacles, but is true for some ( 4, 43. 50 ). For example, iin — 2 and D is 
an ellipse, then the series analogous to H6.6|l converges in the region r > a, where 2a 
is the distance between the foci of the ellipse T . In the engineering literature there 
are numerical algorithms, based on the Rayleigh conjecture. Our aim is to give a 
formulation of a Modified Rayleigh Conjecture (MRC) which holds for any Lipschitz 
obstacle and can be used in numerical solution of the direct and inverse scattering 
problems (see 69 ). We discuss the Dirichlet condition but similar argument is 
applicable to the Neumann boundary condition, corresponding to acoustically hard 
obstacles. 

Fix e > 0, an arbitrary small number. 

Lemma 6.1. There exist L = L{e) and ci = Ci{e) such that 

Lie) 

(6-7) \\uo + Y^ Q (£)'(/'£ 1 1 L2 (5) < e. 

i=0 

If | |6'. ?| ) and the boundary condition 1)6'. J |) hold, then 

Lie) 

(6.8) - ""1^2(5) < e, :=^ce{e)il;e. 

i=o 

Lemma 6.2. // 1)6'. ^) holds then 

(6.9) - will = 0(e), e->0, 

where \\\ ■ \\\ || • \\hi^JD') + II • \\L'^iD'-ii+M)-'), j > 1, m > is an arbitrary 
integer, iJ™ is the Sobolev space, and Ve,v in 1^6. y\} are functions defined in D' . 
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In particular, implies 

(6.10) \\v,-v\\L2iSn.)=0{e), e^O, 
where Sr is the sphere centered at the origin with radius R. 

Lemma 6.3. One has: 

(6.11) Q(e)^A,(a), V^, e^O. 

The Modified Rayleigh Conjecture (MRC) is formulated as a tfieorem, wliich 
follows from the above three lemmas: 

Theorem 6.4. For an arbitrary small e > there exist L(e) and Q(e), < 
e < L{e), such that and Um]) hold. 

See |69| for a proof of the above statements. 

The difference between RC and MRC is: H6.8|l does not hold if one replaces 
by '^i^i)Ai{a)^i, and lets L oo (instead of letting e 0). Indeed, the series 
^^^Q (a)^£ diverges at some points of the boundary for many obstacles. Note 
also that the coefficients in H6.8II depend on e, so H6.8II is not a partial sum of a 
series. 

For the Neumann boundary condition one minimizes 

dN 

L^S) 

with respect to q. Analogs of Lemmas 16. lltOl are valid and their proofs are essen- 
tially the same. 

See [7B for an extension of these results to scattering by periodic structures. 

6.3. Numerical Experiments. In this section we desribe numerical results 
obtained by the Random Multi-point MRC method for 2D and 3D obstacles. We 
also compare the 2D results to the ones obtained by our earlier method introduced in 
|30 |. The method that we used previously can be described as a Multi-point MRC. 
Its difference from the Random Multi-point MRC method is twofold: It is just the 
first iteration of the Random method, and the interior points Xj, j = 1,2,..., J 
were chosen deterministically, by an ad hoc method according to the geometry of 
the obstacle D. The number of points J was limited by the size of the resulting 
numerical minimization problem, so the accuracy of the scattering solution (i.e. 
the residual r™™) could not be made small for many obstacles. The method was 
not capable of treating 3D obstacles. These limitations were removed by using the 
Random Multi-point MRC method. As we mentioned previously, 30 contains a 
favorable comparison of the Multi-point MRC method with the BIEM, inspite in 
spite of the fact that the numerical implementation of the MRC method in |30| is 
considerably less efficient than the one presented in this paper. 

A numerical implementation of the Random Multi-point MRC method follows 
the same outline as for the Multi-point MRC, which was described in [3U |. Of 
course, in a 2D case, instead of H6.4|) one has 

where {x — Xj)/\x — Xj \ — e'^^' . 
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For a numerical implementation choose M nodes {tm} on the surface S of the 
obstacle D. After the interior points Xj^ j ~ 1, 2, J are chosen, form N vectors 

n = 1, 2, . . . , iV of length M. Note that iV = (2L + 1) J for a 2D case, and N = 
(L + 1)^ J for a 3D case. It is convenient to normalize the norm in M*^ by 

1 

m— 1 

Then IImoII = 1- 

Now let b — {g{tm)}m=iJ the Random Multi-point MRC (see section 1), 
and minimize 

(6.12) $(c) = ||b + Ac|l, 

for c e C^, where A is the matrix containing vectors a("), n = 1,2, . . . ,iV as its 
columns. 

We used the Singular Value Decomposition (SVD) method (see e.g. |47) ) to 
minimize ()6.12|) . Small singular values s„ < Wmin of the matrix A are used to 
identify and delete linearly dependent or almost linearly dependent combinations 
of vectors a^"^ . This spectral cut-off makes the minimization process stable, see the 
details in [30| . 

Let r™™ be the residual, i.e. the minimal value o/ $(c) attained after Nmax 
iterations of the Random Multi-point MRC method (or when it is stopped). For a 
comparison, let r^i™ be the residual obtained in |30| by an earlier method. 

We conducted 2D numerical experiments for four obstacles: two ellipses of 
different eccentricity, a kite, and a triangle. The M=720 nodes tm were uniformly 
distributed on the interval [0, 27r], used to parametrize the boundary S. Each case 
was tested for wave numbers k — 1.0 and k = 5.0. Each obstacle was subjected to 
incident waves corresponding to a = (1.0,0.0) and a = (0.0, 1.0). 

The results for the Random Multi-point MRC with J = 1 are shown in Table 
1101 in the last column r™™. In every experiment the target residual e = 0.0001 was 
obtained in under 6000 iterations, in about 2 minutes run time on a 2.8 MHz PC. 

In 3Q|, we conducted numerical experiments for the same four 2D obstacles 
by a Multi-point MRC, as described in the beginning of this section. The interior 
points Xj were chosen differently in each experiment. Their choice is indicated in 
the description of each 2D experiment. The column J shows the number of these 
interior points. Values L = 5 and M = 720 were used in all the experiments. These 
results are shown in Table [TUI column rJJ^™. 

Thus, the Random Multi-point MRC method achieved a significant improve- 
ment over the earlier Multi-point MRC. 

Experiment 2D-I. The boundary S is an ellipse described by 

(6.13) r(f) = (2.0cosi, sini), < t < 27r . 

The Multi-point MRC used J = 4 interior points Xj = 0.7r(^!i^), j = 1, . . . ,4. 
Run time was 2 seconds. 

Experiment 2D-II. The kite-shaped boundary S (see |17| . Section 3.5) is 
described by 

(6.14) r(i) = (-0.65 + cosi + 0.65cos2i, l.Ssini), < t < 2tt . 
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Table 10. Normalized residuals attained in the numerical exper- 
iments for 2D obstacles, ||uo|| = 1. 
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The Multi-point MRC used J = 16 interior points Xj = 0.9r{^^\^), j = 1, . . . , 16. 
Run time was 33 seconds. 

Experiment 2D-III. The boundary S is the triangle with vertices (—1.0, 0.0) 
and (1.0, ±1.0). The Multi-point MRC used the interior points Xj = 0.9r(^^^^^), 
j = 1, . . . , 16. Run time was about 30 seconds. 

Experiment 2D-IV. The boundary S is an ellipse described by 

(6.15) r{t) = (0.1 cost, sini), < t < 27r . 

The Multi-point MRC used J = 32 interior points x^ = 0.95r(^^^ig^), j = 
1, . . . , 32. Run time was about 140 seconds. 

The 3D numerical experiments were conducted for 3 obstacles: a sphere, a cube, 
and an ellipsoid. We used the Random Multi-point MRC with L = 0, Wmin = 
10~^^, and J — 80. The number M of the points on the boundary S is indicated in 
the description of the obstacles. The scattered field for each obstacle was computed 
for two incoming directions ai = {0, 0), z = 1, 2, where </> was the polar angle. The 
first unit vector ai is denoted by (1) in Table [TTl ai — (0.0, 7r/2). The second one is 
denoted by (2), a2 = (7r/2, 7r/4). A typical number of iterations Nuer and the run 
time on a 2.8 MHz PC are also shown in Table [TTl For example, in experiment I 
with k = 5.0 it took about 700 iterations of the Random Multi-point MRC method 
to achieve the target residual r""" = 0.001 in 7 minutes. 

Experiment 3D-I. The boundary S is the sphere of radius 1, with M — 450. 

Experiment 3D-II. The boundary S is the surface of the cube [—1,1]'^ with 
M = 1350. 

Experiment 3D-III. The boundary S is the surface of the ellipsoid x'^/16 -I- 
y2 + = 1 with M = 450. 
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Table 1 1 . Normalized residuals attained in the numerical exper- 
iments for 3D obstacles, ||uo|| = 1. 
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In the last experiment the run time could be reduced by taking a smaller value 
for J. For example, the choice of J = 8 reduced the running time to about 6-10 
minutes. 

Numerical experiments show that the minimization results depend on the choice 
of such parameters as J, Wmin, and L. They also depend on the choice of the in- 
terior points Xj. It is possible that further versions of the MRC could be made 
more efficient by finding a more efficient rule for their placement. Numerical exper- 
iments in j30j showed that the efficiency of the minimization greatly depended on 
the deterministic placement of the interior points, with better results obtained for 
these points placed sufficiently close to the boundary S of the obstacle D, but not 
very close to it. The current choice of a random placement of the interior points 
Xj reduced the variance in the obtained results, and eliminated the need to pro- 
vide a justified algorithm for their placement. The random choice of these points 
distributes them in the entire interior of the obstacle, rather than in a subset of it. 

6.4. Conclusions. For 3D obstacle Rayleigh's hypothesis (conjecture) says 
that the acoustic field u in the exterior of the obstacle D is given by the series 
convergent up to the boundary of D: 

oo 

(6.16) uix,a)=e''"'-^ +Y,Mc^)^e, i^i := Yeia')he{kr), a' = -. 

e=o ^ 

While this conjecture (RC) is false for many obstacles, it has been modified in j69) 
to obtain a valid representation for the solution of H6.1|) - (|6.2I) . This representation 
(Theorem 16. 4|l is called the Modified Rayleigh Conjecture (MRC), and is, in fact, 
not a conjecture, but a Theorem. 

Can one use this approach to obtain solutions to various scattering problems? 
A straightforward numerical implementation of the MRC may fail, but, as we show 
here, it can be efficiently implemented and allows one to obtain accurate numerical 
solutions to obstacle scattering problems. 

The Random Multi-point MRC algorithm was successfully applied to various 
2D and 3D obstacle scattering problems. This algorithm is a significant improve- 
ment over previous MRC implementation described in |30| . The improvement is 
achieved by allowing the required minimizations to be done iteratively, while the 
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previous methods were limited by the problem size constraints. In |30| . such MRC 
method was presented, and it favorably compared to the Boundary Integral Equa- 
tion Method. 

The Random Multi-point MRC has an additional attractive feature, that it can 
easily treat obstacles with complicated geometry (e.g. edges and corners). Unlike 
the BIEM, it is easily modified to treat different obstacle shapes. 

Further research on MRC algorithms is conducted. It is hoped that the MRC 
in its various implementation can emerge as a valuable and efficient alternative to 
more established methods. 

7. Support Function Method for inverse obstacle scattering problems. 

7.1. Support Function Method (SFM). The Inverse Scattering Problem 
consists of finding the obstacle D from the Scattering Amplitude, or similarly ob- 
served data. The Support Function Method (SFM) was originally developed in a 
3-D setting in j48| . see also |50j . pp 94-99. It is used to approximately locate the 
obstacle D. The method is derived using a high-frequency approximation to the 
scattered field for smooth, strictly convex obstacles. It turns out that this inexpen- 
sive method also provides a good localization of obstacles in the resonance region 
of frequencies. If the obstacle is not convex, then the SFM yields its convex hull. 

One can restate the SFM in a 2-D setting as follows (see |31p. Let D C be 
a smooth and strictly convex obstacle with the boundary F. Let ^'(y) be the unique 
outward unit normal vector to F at y G F. Fix an incident direction a G S^. Then 
the boundary F can be decomposed into the following two parts: 

(7.1) F+ = {y e F : u{y) • a < 0} , and F_ = {y e F : iy(y) • a > 0} , 

which are, correspondingly, the illuminated and the shadowed parts of the boundary 
for the chosen incident direction a. 

Given a G S^, its specular point So(q;) G F+ is defined from the condition: 

(7.2) So(a) • a — min s • a 

Note that the equation of the tangent line to F+ at Sq is 

(7.3) < xi,X2 > ■ a = So{a) ■ a , 
and 

(7.4) z/(so(a)) = -a. 
The Support function d{a) is defined by 

(7.5) d{a) — So(a) • a . 

Thus |d(Q;)| is the distance from the origin to the unique tangent line to F+ 
perpendicular to the incident vector a. Since the obstacle D is assumed to be 
convex 

(7.6) 1) = n„gsi{x e : x-a>d(a)}. 

The boundary F of Z? is smooth, hence so is the Support Function. The knowl- 
edge of this function allows one to reconstruct the boundary F using the following 
procedure. 

Parametrize unit vectors 1 € 5^ by l{t) = (cost, sin i), < t < 2tt and define 

(7.7) p{t) ^ d(l{t)), 0<t<27r. 
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Equation (|7.,S|I and the definition of the Support Function give 

(7.8) xi cos t + X2sait = p{t) . 

Since F is the envelope of its tangent hues, its equation can be found from (|7.8() 
and 

(7.9) —xi sin i + a;2 cost = p'(<) . 
Therefore the parametric equations of the boundary F are 

(7.10) xi[t) — p{t) cost ~ p' {t)su\t, X2{t) — p{t) sint + p'{t) cost . 

So, the question is how to construct the Support function d{V), I E from the 
knowledge of the Scattering Amphtude. In 2-D the Scattering Amphtude is related 
to the total field u = uq + v hy 

(7.11) A{a', a) ^ [ e-^-''^ d.(y) . 

In the case of the "soft" boundary condition (i.e. the pressure field satisfies the 
Dirichlet boundary condition u — 0) the Kirchhoff (high frequency) approximation 
gives 

(7.12) '^^2^ 

on the illuminated part F_|_ of the boundary F, and 

(7.13) 1^ = 

on the shadowed part F_. Therefore, in this approximation, 

(7.14) A(a',a) = -^^ / a • i/(y) e*'=("-"'>y ds(y) . 

V 27rfc , 



+ 



Let L be the length of F+, and y = y(C)j < C, 1^ L he \is arc length 
parametrization. Then 



L 



(7.15) A{a',a) = / a ■ v{y{0) eifc("-"')-y(C) . 

\/27r 



Let C,Q e [0, L] be such that Sq = y(Co) is the specular point of the unit vector 
I, where 



(7.16) 1 



\a — a'\ 

Then v{so) = -1, and d(\) = y(Co) • 1. Let 

^(0 = (a - a') ■ y(C) . 
Then tf{C,) — 1 • y(C)|a — a'\. Since j^(so) and y'(Co) are orthogonal, one has 

^'(Co) = i-y'(Co)l«-a'l = o. 

Therefore, due to the strict convexity of D, C,q is also the unique non-degenerate 
stationary point of <p(C) on the interval [0, L], that is ^p' {C,o) = 0, and <p"(Co) 7^ 0. 
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According to the Stationary Phase method 
(7.17) 



/(C)' 



dC = /(Co) exp 



(^"(Co) 
4 \v"iCo)\ 



2tt 



k\v"{Co)\ 



l + O 



as /c ^ 00. 

By the definition of the curvature k(Co) ~ |y"(Co)|- Therefore, from the 
colhnearity of y"(Co) and 1, |<y3"(Co)| — \a — q;'|k(Co)- Finally, the strict convex- 
ity of D, and the definition of (y9(C), imply that Co is the unique point of minimum 
of on [0, L], and 

Using H7.17|l - H7.18|l . expression (|7.15|) becomes: 



(7.19) v4(a',a) 



1-a 



Jk{a-a')-y{Co) 



1 + 



y^\a - a'|K(Co) 

At the specular point one has l-a' = — l-a. By the definition a — a' — Ija — q;'|. 
Hence 1- {a — a') = ja — q;'| and 21- a — \a — a'\. These equalities and d(l) — y(Co) '1 
give 



(7.20) A{a',a)^-^^ 
Thus, the approximation 
(7.21) 



«(Co) 



Jk\a-a'\d{l} 



l + O 



oo . 



A{a',a) « --^ 



^ik\a~a.' \d(l) 



can be used for an approximate recovery of the curvature and the support function 
(modulo 27r/fc|a — a'\) of the obstacle, provided one knows that the total field 
satisfies the Dirichlet boundary condition. The uncertainty in the support function 
determination can be remedied by using different combinations of vectors a and a' 
as described in the numerical results section. 

Since it is also of interest to localize the obstacle in the case when the boundary 
condition is not a priori known, one can modify the SFM as shown in [T7| . and 
obtain 



(7.22) 
where 

and 



A{a',a) 



«(Co) 



„i(fc|Q-a'|d(l)-27o+Tr) 



7o = arctan — , 
h 



9u 

— + hu = 
on 



along the boundary T of the sought obstacle. 

Now one can recover the Support Function d{l) from H7.22II . and the location 
of the obstacle. 
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7.2. Numerical results for the Support Function Method. In the first 
numerical experiment the obstacle is the circle 



(7.23) 



D ^ {{xi,X2) eM.'^ : (xi - 6)2 + (a;2 - 2)2 = 1} , 



It is reconstructed using the Support Function Method for two frequencies in the 
resonance region: k = 1.0, and k = 5.0. Table ^] shows how well the approxi- 
mation H7.21|l is satisfied for various pairs of vectors a and a' all representing the 
same vector 1 = (1.0,0.0) according to H7.16|l . The Table shows the ratios of the 
approximate Scattering Amplitude Aa{a',a) defined as the right hand side of the 
equation (|7.21|l to the exact Scattering Amplitude A{a' ,a). Note, that for a sphere 
of radius a, centered at Xq GR^, one has 



(7.24) 



A{a', a) = 



nk ^ 



Ji{ka) 



H['\ka) 



where a' ~ x/|x| — e'^, and a = e^^. Vectors a and a' are defined by their polar 
angles shown in Table IT^ 

Table 12. Ratios of the approximate and the exact Scattering 
Amplitudes Aa{a' ,a)/A(a' ,a) for 1 = (1.0,0.0). 



a' 


a 


k 


= 1.0 


k 


= 5.0 


TV 





0.88473 - 


0.17487i 


0.98859 - 


0.05846i 


237r/24 


7r/24 


0.88272 - 


0.17696i 


0.98739 - 


0.06006i 


227r/24 


27r/24 


0.87602 - 


0.18422i 


0.98446 - 


0.06459i 


2l7r/24 


37r/24 


0.86182 - 


0.19927j 


0.97977 - 


0.07432i 


207r/24 


47r/24 


0.83290 - 


0.22411j 


0.96701 - 


0.08873i 


197r/24 


57r/24 


0.77723 - 


0.25410j 


0.95311 - 


0.10321i 


187r/24 


67r/24 


0.68675 - 


0.27130j 


0.92330 - 


0.14195i 


177r/24 


77r/24 


0.57311 - 


0.25360j 


0.86457 - 


0.14959i 


167r/24 


87r/24 


0.46201 - 


0.19894j 


0.81794 - 


0.22900i 


157r/24 


97r/24 


0.36677 - 


0.12600j 


0.61444 - 


0.19014i 


147r/24 


107r/24 


0.28169 - 


0.05449j 


0.57681 - 


0.31075i 


137r/24 


ll7r/24 


0.19019 + 


0.00075j 


0.14989 - 


0.09479i 


127r/24 


127r/24 


0.00000 + 


O.OOOOOi 


0.00000 + 


O.OOOOOi 



Table El shows that only vectors a close to the vector 1 are suitable for the 
Scattering Amplitude approximation. This shows the practical importance of the 
backscattering data. Any single combination of vectors a and a' representing 1 is 
not sufficient to uniquely determine the Support Function (i(l) from H7.21|l because 
of the phase uncertainty. However, one can remedy this by using more than one 
pair of vectors a and a' as follows. 

Let 1 e 5^ be fixed. Let 



R{1) = {a e 5^ : |a • 1| > 1/V2}. 
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Define * 



R+ by 



A{a',a) 



-a'\t 



\Aia',a)\ 



where a' — a' (a) is defined by 1 and a according to (|7.1t)|) . and the integration is 
done over a G 

If the approximation (|7.21() were exact for any a € -R(l), then the value of 
^{d(l)) would be zero. This justifies the use of the minimizer to G R of the function 
^'(t) as an approximate value of the Support Function d(l). If the Support Function 
is known for sufficiently many directions 1 S 5*^ , the obstacle can be localized using 
(17. 61) or H7.10|l . The results of such a localization for fc = 1.0 together with the 
original obstacle D is shown on Figure[S| For k = 5.0 the identified obstacle is not 
shown, since it is practically the same as D. The only a priori assumption on D 
was that it was located inside the circle of radius 20 with the center in the origin. 
The Support Function was computed for 16 uniformly distributed in vectors 1. 
The program run takes about 80 seconds on a 333 MHz PC. 



4 



2 




2 4 6 8 

Figure 5. Identified (dotted line), and the original (solid line) 
obstacle D for k = 1.0. 

In another numerical experiment we used A: = 1.0 and a kite-shaped obstacle. 
Its boundary is described by 

(7.25) r(i) = (5.35 + cost + 0.65cos2i, 2.0 + 1. 5sint), < i < 27r . 

Numerical experiments using the boundary integral equation method (BIEM) for 
the direct scattering problem for this obstacle centered in the origin are described 
in |17| . section 3.5. Again, the Dirichlet boundary conditions were assumed. We 
computed the scattering amplitude for 120 directions a using the MRC method 
with about 25% performance improvement over the BIEM, see |30| . 



44 



ALEXANDER G. RAMM AND SEMION GUTMAN 



y 



4 



2 




X 



2 



4 



6 



8 



Figure 6. Identified points and the original obstacle D (solid 
line); k = 1.0. 

The Support Function Method (SFM) was used to identify the obstacle D 
from the synthetic scattering amplitude with no noise added. The only a priori 
assumption on D was that it was located inside the circle of radius 20 with the center 
in the origin. The Support Function was computed for 40 uniformly distributed in 
vectors 1 in about 10 seconds on a 333 MHz PC. The results of the identification 
are shown in Figure El The original obstacle is the solid line. The points were 
identified according to H7.10|l . As expected, the method recovers the convex part of 
the boundary F, and fails for the concave part. The same experiment but with k — 
5.0 achieves a perfect identification of the convex part of the boundary. In each case 
the convex part of the obstacle was successfully localized. Further improvements in 
the obstacle localization using the MRC method are suggested in ,69) . and in the 
next section. 

For the identification of obstacles with unknown boundary conditions let 



where, given t, the vectors a and a' are chosen as above, and the phase function 
il){t), \pl < i < 2 is continuous. Similarly, let v4a(i), i/'al^) be the approximate 
scattering amplitude and its phase defined by formula H7.22|l . 

If the approximation H7.22|l were exact for any a e i?(l), then the value of 



would be a multiple of Itx. 

This justifies the following algorithm for the determination of the Support 
Function (i(l): 

Use a linear regression to find the approximation 



|V'a(i)-fcMl) + 27o-7r| 



on the interval \f2 < t < 2. Then 



(7.26) 
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Table 13. Identified values of the Support Function for tlie circle 
of radius 1.0 a.t k — 3.0. 



h 


Identified d{l) 


Actual d{l) 


0.01 


-0.9006 


-1.00 


0.10 


-0.9191 


-1.00 


0.50 


-1.0072 


-1.00 


1.00 


-1.0730 


-1.00 


2.00 


-0.9305 


-1.00 


5.00 


-1.3479 


-1.00 


10.00 


-1.1693 


-1.00 


100.00 


-1.0801 


-1.00 



Also 

h = — fc tan — - . 

2 

However, the formula for h did not work well numerically. It could only determine 
if the boundary conditions were or were not of the Dirichlet type. Table [T51 shows 
that the algorithm based on 1)7. 26|) was successful in the identification of the circle of 
radius 1.0 centered in the origin for various values of h with no a priori assumptions 
on the boundary conditions. For this circle the Support Function d{l) — —1.0 for 
any direction 1. 

8. Analysis of a Linear Sampling method. 

During the last decade many papers were published, in which the obstacle iden- 
tification methods were based on a numerical verification of the inclusion of some 
function / := /(a, z), z G M'^, a G S'^, in the range R{B) of a certain operator B. 
Examples of such methods include |14| . jl6 | .|42 | . However, one can show that the 
methods proposed in the above papers have essential difficulties, see '76'. Although 
it is true that / ^ R{B) when z ^ D, it turns out that in any neighborhood of / 
there are elements from R{B). Also, although / e R{B) when z E D, there are 
elements in every neighborhood of / which do not belong to R{B) even if z G D. 
Therefore it is quite difficult to construct a stable numerical method for the iden- 
tification of D based on the verification of the inclusions / ^ R{B), and / € R{B). 
Some published numerical results were intended to show that the method based 
on the above idea works practically, but it is not clear how these conclusions were 
obtained. 

Let us introduce some notations : N(B) and R{B) are, respectively, the null- 
space and the range of a linear operator B, D EM.^^ is a. bounded domain (obstacle) 
with a smooth boundary S, D' = \ D, uq — e^^""'^ , k = const > 0, a G S*^ 
is a unit vector, N is the unit normal to S pointing into D' , g — g(x,y,k) := 
g{\x - y\) := , / := e"'''"'''', where z G and a' G 5^, a' := xr~^, r = |a;|, 

u ~ u(x, a, k) is the scattering solution: 

(8.1) (A + fc2)M = in D',u\s^O, 

(8.2) u~uq + v, V — A{a\a,k)e''^'^r^^ + o{r^^), as r ^ oo. 
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where A := A{a' , a, k) is called the scattering amplitude, corresponding to the ob- 
stacle D and the Dirichlet boundary condition. Let G — G{x, y, k) be the resolvent 
kernel of the Dirichlet Laplacian in D': 

(8.3) {A + k^)G^-5{x-y) in D',G|5 = 0, 

and G satisfies the outgoing radiation condition. 
If 

(8.4) (A + fc2)w = in D',w\s = h, 
and w satisfies the radiation condition, then ( 50 ) one has 

(8.5) w{x)^ ( GN{x,s)h{s)ds, w ^ A{a' ,k)e'''''r-^ + o{r-^), 

Js 

as r — > cxD, and xr~^ — a' . We write A{a') for A{a' , k), and 

(8.6) A{a') -.^ Bh -.^ — [ UN(s,-a')h(s)ds, 

47r Js 

as follows from Ramm's lemma: 

Lemma 1. (|50). p. 46) One has: 

(8.7) G{x,y,k)=g{r)u{y,—a\k)+o{r~^), as r=\x\^oo, xr~^ — a' , 

where u is the scattering solution of f^.lp - f^.l^} . 
One can write the scattering amplitude as: 

(8.8) Aia',a,k) = [ ujv(s, -a')e'''"''ds. 

47r Js 

The following claim follows easily from the results in |5Uj . |55) (cf |42) ) : 
Claim: f := e"**^"'-^ £ R{B) if and only if z G D. 
Proof: If e"'''" ~ Eh, then Lemma 1 and (12.6) imply 

9{v^z)= / GN{s,y)hds for ly] > |z| . 
Js 

Thus z Cz D, because otherwise one gets a contradiction: lim^^^ g{y, z) = oo if 
z Cz D' , while lim^^^ Jg Gn{s, y)hds < oo if z € D' . Conversely, if z S Z), then 
Green's formula yields g{y,z) = GN{s,y)g{s, z)ds. Taking \y\ oo, = a', 

and using Lemma 1, one gets e^'*^" = Bh, where h — g{s, z). The claim is proved. 
□ 

Consider B : L^{S) L^{S^), and A : L^{S^) L^{S^), where B is defined 
in (|8.6|l and Aq :— /^a A{a' ,a)q{a)da. Then one proves (see |76j): 
Theorem 1. The ranges R{B) and R{A) are dense in L^{S^) 
Remark 1. In 14 the 2D inverse obstacle scattering problem is considered. 
It is proposed to solve the equation (1.9) in |14| : 

(8.9) / A{a,f3)gdf3 ^ e-''"^-', 

Js^ 

where A is the scattering amplitude at a fixed k > 0, is the unit circle, a G S^, 
and z is a point on R^. If ^ = Q{(3,z) is found, the boundary S of the obstacle 
is to be found by finding those z for which \\Q\\ :— \\Q{(3, z)\\l2(^s^) is maximal. 
Assuming that k'^ is not a Dirichlet or Neumann eigenvalue of the Laplacian in D, 
that D is a smooth, bounded, simply connected domain, the authors state Theorem 
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2.1 jl4j . p. 386, which says that for every e > there exists a function Q e L^{S^), 
such that 

(8.10) Urn ||e(/3,z)|| =oo, 

z—>S 

and (see [H], p.386), 

(8.11) II / A{a,p)gdP-e-'''^-'\\<e. 

There are several questions concerning the proposed method. 

First, equation (|8.9|l . in general, is not solvable. The authors propose to solve 
it approximately, by a regularization method. The regularization method applies 
for stable solution of solvable ill-posed equations (with exact or noisy data). If 
equation H8.9|l is not solvable, it is not clear what numerical "solution" one seeks 
by a regularization method. 

Secondly, since the kernel of the integral operator in ((8.9(1 is smooth, one can 
always find, for any z G R^, infinitely many Q with arbitrary large \ \G\\, such that 
((8.11(1 holds. Therefore it is not clear how and why, using ((8.10(1 . one can find S 
numerically by the proposed method. 

A numerical implementation of the Linear Sampling Method (LSM) suggested 
in [14[ consists of solving a discretized version of 1(8.9(1 

(8.12) Fg - f , 

where F = {Aai,Pj}, i = 1,...,N, j = 1,...,N be a square matrix formed by 
the measurements of the scattering amplitude for N incoming, and N outgoing 
directions. In 2-D the vector f is formed by 



see |10| for details. 

Denote the Singular Value Decomposition of the far field operator by = 
USV^. Let Sn be the singular values of F, p = U"i, and ^ = V"{. Then the 
norm of the sought function g is given by 

(8.13) ||t;f^^^. 

n=l " 

A different LSM is suggested by A. Kirsch in ,42j. In it one solves 

(8.14) {F*Fy/^g = { 

instead of 1(8.12(1 . The corresponding expression for the norm of Q is 

N I |2 

(8.15) ||t;|p = ^^. 

71— 1 

A detailed numerical comparison of the two LSMs and the linearized tomographic 
inverse scattering is given in [1U[ . 

The conclusions of 'lUI, as well as of our own numerical experiments are that 
the method of Kirsch 1(8.14(1 gives a better, but a comparable identification, than 
(18.12(1 . The identification is significantly deteriorating if the scattering amplitude is 
available only for a limited aperture, or the data are corrupted by noise. Also, the 
points with the smallest values of the \\G\\ are the best in locating the inclusion, and 
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gnk, k.=1.0, noise=0.00 




Figure 7. Identification of a circle at fc = 1.0. 



not tfie largest one, as required by the theory in |42| and in |14| . In Figures [7| and |S| 
the implementation of the Colton-Kirsch LSM (|8.13|l is denoted by gnck, and of the 
Kirsch method H8.15|l by gnk. The Figures show a contour plot of the logarithm of 
the II 5 II . In all the cases the original obstacle was the circle of radius 1.0 centered 
at the point (10.0, 15.0). A similar circular obstacle that was identified by the 
Support Function Method (SFM) is discussed in Section 10. Note that the actual 
radius of the circle is 1.0, but it cannot be seen from the LSM identification. The 
LSM does not require any knowledge of the boundary conditions on the obstacle. 
The use of the SFM for unknown boundary conditions is discussed in the previous 
section. The LSM identification was performed for the scattering amplitude of the 
circle computed analytically with no noise added. In all the experiments the value 
for the parameter TV was chosen to be 128. 
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